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Abstract and Keywords 

Abstract: This dissertation explored topographic controls on spatial and temporal 
patterns in water yield and nutrient (carbon, nitrogen and phosphorus) export from 
forested headwater catchments in the Turkey Lakes Watershed in central Ontario, where 
other factors contributing to differences in water yield and nutrient export, including 
climate, geology, forest, and soils, are relatively constant. Topographic characteristics, 
including (a) hydrological flushing potential (expansion of water table into nitrate-N 
producing areas); (b) hydrological storage potential (area of wetlands, which can 
alternatively allow water and nutrients to bypass wetlands when storage capacity is filled 
with water or to trap them when not filled); and (c) hydrological loading potential 
(differences in precipitation caused by elevation), were considered in deconstructing non- 
stationary (linear trends) and stationary (oscillating cycles) patterns in water yield and 
nutrient export data. Topography explained the majority of differences in water yield and 
nutrient export. For spatial variation, topographic metrics representing hydrologic 
flushing potential predicted the majority of the spatial variation in nitrate-N export. In 
contrast, topographic metrics representing hydrologic storage potential explained the 
majority of the observed spatial variation in dissolved organic carbon, dissolved organic 
nitrogen and total dissolved phosphorus export. For temporal variation, catchments with 
low hydrologic loading potential were generally more sensitive to trends and cycles for 
water and nutrient export. Among these catchments, hydrological storage potential had 
no significant effect on water export trends, but had a significant effect on water export 
cycles; namely, the water export range was larger in the catchments with higher 
hydrological storage potential, even though the water export average was the same as 
catchments with lower hydrological storage potential. For nutrient export, the non- 
stationary signals were not consistent among the nutrients, but the amplitude of stationary 
signals in nutrient export in catchments with high hydrological storage potential 
compared to those with low hydrological storage potential was higher for organic 
nutrients and lower for nitrate-nitrogen. Despite many similarities in these headwater 
catchments, topography influenced the absolute and relative magnitude of hydrological 
and biogeochemical export from these catchments, which will have implications on the 
productivity and biodiversity of downstream aquatic systems. 
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Chapter 1: Introduction 


1.1 Problem Statement 

Despite continued efforts to understand the major controls of nutrient transport on a 
landscape, complete understanding of such controls remains challenging (Neff et al . , 
2003; Deng, et al., 2011) given the need to deal with and characterize spatially 
heterogeneous landscapes and temporally dynamic climate related transport controls 
(Stieglitz et al., 2003; Christopher et al., 2008). Heterogeneity in landscape 
characteristics not only makes identification of potential sources and sinks of nutrients on 
a landscape difficult but also results in complexity in understanding the net effect of 
interactions between potential sources and sinks of nutrients. Furthermore, climate is 
temporally dynamic and often contains a complex mixture of non- stationary linear trends 
and nested stationary signals of a cyclical nature, making the investigation of climatic 
control on nutrient transport challenging. Therefore, there is need for improved 
characterization and understanding of spatial and temporal drivers of nutrient transport on 
a landscape, because such understanding is central to implementing effective water 
quality protection strategies (Neal and Heathwaite, 2005; Robinson et al., 2005). 

1.2 Scientific Rationale 

Many hydrological and biogeochemical processes can only be understood at hillslope or 
small experimental catchment scales (Harris and Heathwaite, 2005; McDonnell et al., 
2007). In the last few decades, catchments have been used as important units of 
investigation in the study of biogeochemical processes and examination of the spatial and 
temporal controls in nutrient transport. A catchment is an area of land that funnels water 
and associated sediment and nutrients to the lowest point within its boundary, the 
catchment outlet. Here, water as well as the associated nutrients and sediment are 
exported outside of the catchment’s boundary to receiving water bodies, integrating the 
effects of nutrient transport controls such as climate, landscape characteristics and land- 
use changes (Ocampo et al., 2006a,b; Saunders et al., 2006; Creed and Beall, 2009; 

Piatek et al., 2009; Mitchell, 201 1). Particularly, headwater catchments, being small in 
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size and having thin soil layers,, steep slopes and low storage, can quickly respond to 
small environmental changes and easily detect spatial differences in landscape 
characteristics, and therefore have the potential to serve as sentinels for climate 
associated temporal changes and/or landscape related spatial heterogeneities. Strahler 
(1957) classified a headwater catchment as an area of land draining to a first or second- 
order stream channel, even though, such classification is controversial because the 
resolution of the scale of map used for the definition of stream channels affects the 
identification and interpretation of the size of headwater catchments (Gomi et al., 2002). 

The deterioration of water quality due to increased nutrient deposition results in 
high biological productivity (eutrophication), which has numerous environmental, social 
and economic implications (Prepas et al., 1997; Kangur and Mols, 2008; Simis et al., 
2005; Randolph et al., 2008; Dodds et al., 2009). Generally, excessive enrichment of 
surface water bodies with sediment and nutrients (mainly phosphorus and nitrogen) from 
agricultural and/or urban runoff has been the most common explanation for water quality 
deterioration (Karr and Dudley, 1981; Carpenter et al., 1998; Sharpley et al., 2000). 

Recent studies have observed increased deterioration of surface waters in relatively 
pristine forested landscapes such as the Algoma Highlands of Central Ontario (Winter et 
al., 201 1). Mostly, the drivers of biogeochemical dynamics and nutrient movement in 
such pristine forested landscapes are complex and mainly linked to landscape-related 
spatial and climate-associated temporal controls (Williams et al., 1996; Mattsson et al., 
2003; Fisher et al., 2004; Harris and Watmough et al., 2004; Heathwaite, 2005; Creed 
and Beall, 2009). Located in the Canadian Precambrian Shield, the Algoma Highland 
region is dominated by non-carbonate bedrock geology with very poor capacity to buffer 
sulfur and nitrogen deposition in acid rain (Scheider et al., 1979; Shrivastava, 2007). 

Even though sulfur deposition has declined over the last few years, nitrogen continues to 
be the predominant acidifying pollutant in the region (Jeffries et al., 1988; Morrison et al., 
1992; Jeffries, 1995; Jeffries et al., 2003; Creed and Beall, 2009). As the current trend of 
north-ward expansion of human development activities intensifies, the surface water 
quality and consequently the ecosystem integrity of water bodies in such a sensitive 
region can, therefore, be easily threatened without taking appropriate measures that are 
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based on scientific understanding of the controls of biogeochemical dynamics and 
nutrient transport in the region. 

This thesis explores the influence of landscape linked spatial and climate related 
temporal characteristics on hydrological and biogeochemical responses, specifically 
water yield and nitrate-nitrogen (NO 3 -N), dissolved organic nitrogen (DON), dissolved 
organic carbon (DOC) and total dissolved phosphorus (TDP), from small headwater 
catchments in the Turkey Lakes Watershed (TLW) of the Algoma Highlands. These 
nutrients were chosen due to the important role they play in primary productivity. The 
spatial controls on water and nutrient export from headwater catchments include physical 
characteristics such as topography, geology, soil, vegetation and other characteristics that 
are spatially heterogeneous. Heterogeneity in catchment physical characteristics can 
translate to spatial differences in nutrient availability and/or mobility to receiving water 
bodies (Heathwaite and Johnes, 1996; Geyer and Schumann, 2001; Lu et al., 2007; 
Fraterrigo and Downing, 2008; Creed and Beall, 2009), causing catchments to exhibit 
differential hydrologic and biogeochemical responses and reveal differential sensitivity to 
signals of climatic alteration. 

Climate has also an effect on regional hydrology, which in turn influences 
nutrient production and transport, and thus the water quality of receiving surface water 
bodies (Whitehead et al., 2009). Catchment responses to climate are temporally dynamic 
and contain an often complex mixture of non- stationary and nested stationary signals. 
Non-stationary catchment response signals have statistical means and standard deviations 
that change over time. They are mainly linked to unidirectional climate warming driven 
by anthropogenic climate change. Stationary catchment response signals have statistical 
means and standard deviations that do not change over time. They are linked to naturally 
varying large-scale climatic oscillations, caused by fluctuations either in sea surface 
temperature (SST) or sea level pressure (SLP) in areas that can be thousands of 
kilometres away from the site of investigation. By influencing regional climate, climatic 
oscillations can affect hydrological and biogeochemical processes, which in turn may 
have water quality and quantity implications (Keener et al., 2010; Kahya, 2011; 
Niedzielski, 2011). In regions with unidirectional climate change, the effect of climate on 
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hydrological and biogeochemical processes can be significantly manifested when the 
stationary components of climate contribute to the positive phase of oscillations. 

Previous studies investigating climate control on catchment responses did not 
discriminate non- stationary and stationary signals of catchment responses. The majority 
of these studies tried to identify connections between climatic oscillations and stream 
hydrologic responses (e.g., Bonsai and Shabbar, 2008; Burn, 2008; Whitfield et al, 2010; 
Assam et al., 2010), but very little work has focused on relating climatic oscillations to 
catchment nutrient export (Sang et al., 201 1). A scientifically sound, analytical approach 
is needed to discriminate the non- stationary and stationary components of catchment 
hydrological and biogeochemical responses and decompose the nested stationary 
response components into signals of various frequencies of interest. The decomposition 
helps to reveal catchment signals that are masked by the most dominant signal in 
catchment responses. In addition, the decomposed signals can be independently tested 
with climatic oscillations of interest. 

The TLW has long-term water yield and nutrient export data. The catchments are 
relatively pristine, preserving the hydrologically and biogeochemically complex 
landscapes that are necessary to determine a baseline understanding of how the natural 
systems function, while at the same time allowing comparisons with altered systems 
(Mattsson et al., 2003). The watershed is located in a region that shows noticeable annual 
climatic variability and exhibits substantial landscape spatial heterogeneity in terms of 
topography over a small watershed area of only 10.5 km”. In addition, the watershed is 
close to a meteorological station with long-term climatic records. Analysis of the long 
term climatic record in the region indicates that the TLW has been experiencing a 
changing climate, with almost a 1°C per decade rise in average temperature and 80 mm 
per decade decrease in annual precipitation over the past three decades. 

1.3 Goals 

The goal of this study was to improve our understanding of the temporal and spatial 
controls on catchment water yield and nutrient export characteristics by (1) presenting 
simple topographic metrics that represent landscape nutrient sources and sinks, which can 
be used to track nutrient transport and export to stream networks from catchments in 
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forested landscapes, and (2) developing a novel analytical framework for decomposing 
catchment water and nutrient export signals, so that potential for differential sensitivity of 
catchment responses in yearly water and nutrient export to climatic forcing can be 
established. 

1.4 Hypotheses and Questions 

This thesis tests the following hypotheses: (1) topography determines the dominant 
hydrologic processes that influence nutrient transport on a landscape; (2) catchment 
yearly water yield and nutrient export responses contain a complex mixture of non- 
stationary and stationary climate signals, the significance of which is different among 
catchments with different topographic characteristics; (3) non- stationary trends contribute 
high rates over time in nutrient export and can obscure the contribution of stationary 
cycles to nutrient export signal; (4) metabolically more active nitrogen and phosphorus 
containing nutrients reveal stronger decline in their rates of change compared to 
metabolically less active carbon containing nutrients; and (5) catchments with higher 
hydrologic storage and lower hydrologic loading potentials are sensitive to non-stationary 
climatic trends in their yearly nutrient export, mainly due to the influence these 
catchment characteristics have on increasing water residence time and therefore nutrient 
reactivity time. To test these hypotheses, the following three questions were posed: 

Question 1: What relationships are there between spatial topographic metrics and 
nutrient export? 

The study investigated the role of the topographic variation prevalent in the TLW region 
in terms of explaining the observed spatial variation in nutrient export among TLW 
headwater catchments. 

Question 2: Are there differences in the sensitivity of catchment water export 
behaviour to non-stationary climate change and naturally occurring stationary 
climate oscillations? 

To answer this question, four representative headwater catchments were selected based 
on a gradient in water loading (precipitation) and water storage in wetlands. An analytical 
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framework was first developed, and then 28-year time series of water yield were analysed 
for non- stationary (linear trends) and nested stationary (non-linear oscillations) signals. 

Question 3: How are climate forcing and non-stationary and stationary catchment 
nutrient export signals linked? 

The analytical framework used to reveal non-stationary and stationary signals in water 
yields was used to detect signals in yearly time series of nutrient export [dissolved 
organic carbon (DOC), dissolved organic nitrogen (DON), nitrate (NO 3 -N), and total 
dissolved phosphorus (TDP)]. 

1.5 Thesis Organization 

This thesis has been prepared in the integrated article format and consists of three 
manuscripts related to the three main questions. The introduction (Chapter 1) provides an 
overview of the research problem, including the question and objectives that form the 
basis of the thesis. The first manuscript (Chapter 2) investigates spatial controls on 
catchment carbon, nitrogen and phosphorus export patterns. The second manuscript 
(Chapter 3) examines the differential sensitivity of four representative TLW catchments 
to non-stationary and stationary climate signals in their yearly water export responses. A 
detailed explanation of the analytical framework developed for the analysis is provided in 
this chapter. The third manuscript (Chapter 4) applies the analytical framework to reveal 
complex non-stationary and stationary biogeochemical responses in the export of carbon, 
nitrogen, and phosphorus nutrients, as well as investigating the sensitivity of responses to 
climatic forcing and change. The last chapter (Chapter 5) summarizes the major 
conclusions of the work, draws the linkages among the individual chapters, and identifies 
future research directions. 
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Chapter 2: Searching for similarity in topographic 
controls on carbon, nitrogen and phosphorus export 
from forested headwater catchments 


Preface 

Topography influences hydrological processes that in turn affect biogeochemical export 
to surface waters on forested landscapes. In this chapter, the difference in long-term 
average annual dissolved organic carbon (DOC), organic and inorganic nitrogen (NOf- 
N, DON), and phosphorus (TDP) export from catchments with similar climate, geology, 
and forest type and age are established. This chapter tests the hypothesis that topography 
determines the dominance in hydrologic processes that in turn determines nutrient export. 
Topographic indicators are designed to represent topographically-regulated 
hydrological processes that influence nutrient export: (a) hydrological storage potential 
( i. e. , effect of topographic flats/depressions in a catchment, representing wetlands that 
may either result in increased runoff when the water table in the wetland is near or at the 
surface, or decreased runoff when the water table in the wetland is well below the surface 
preventing water contributing to the wetland from being exported to the stream); and (b) 
hydrological flushing potential (i.e., effect of topographic slopes on potential for variable 
source area to expand and tap into previously untapped nutrient source areas in the 
forest floor and surface mineral layers). These topographic indicators are then related to 
nutrient export, to see if shifting dominance in hydrological processes may explain 
differences in nutrient export. The significance of this chapter is that it presents simple 
topographic indicators that can be used to track nutrient sources, sinks, and their 
transport and export to surface waters from catchments on forest landscapes. 

2.1 Introduction 

General rules of nutrient export derived from careful study of gauged catchments are 
required to make informed predictions about responses in ungauged basins. This is the 
first step towards calculating cumulative impacts of terrestrial-atmospheric-aquatic 
linkages in larger basins. Within landscapes where climate, geology, and biology are 
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similar, topography is the key control on nutrient accumulation, transformation and 
transport through its fundamental influence on energy and water mass balances (Sommer 
and Schlichting, 1997). Topographic characteristics, which are measured or 
mathematically analysed configurations of the landform shape and structure, influence 
the flow of water on the landscape (Wagener et al., 2007). In turn, hydrologic flow paths 
influence the transport and transformation of nutrients. Several studies have shown strong 
associations between topographic properties and nutrient export behavior of carbon 
(Creed et al., 2003, 2008), nitrogen (Creed and Beall, 2009) and phosphorus (Devito el 
al., 2000). However, the nature of topographic controls on nutrient dynamics is complex, 
suggesting that the dominant hydrologic process influencing export may be different for 
each nutrient. 

Topography is important for understanding the aquatic fate of dissolved carbon, 
nitrogen and phosphorus. For example, following a rain event on landscapes where 
topography is controlled by bedrock, the groundwater table rises toward the surface and 
intersects with surface soils that have accumulated nutrients in the intervening dry period. 
These nutrients are then mobilized and flushed to receiving waters resulting in export of 
carbon (Homberger et al., 1994), nitrogen (Creed et al., 1996) or phosphorus (Evans et 
al., 2000). In these cases topography influences the hydrological flushing through 
transport of nutrients to surface waters as a function of the size and spatial organization 
of the variable source area. Topography provides a conduit for nutrient export, initially 
through slower shallow, subsurface flows and then more rapidly as nutrients reach zones 
of saturation directly connected to streams. However, there are instances in which the 
topography may result in nutrients being transformed. For example, wetlands, 
depressions or flat areas intercept nutrient runoff. The result of this interception may 
either reduce the amount of nutrient export from catchments (e.g., NOf-N. Creed et al., 
2009) or increase the supply of nutrients for export [e.g., carbon, phosphorus (Creed et al., 
2008) and organic nitrogen (Creed and Beall, 2009)]. 

Linking hydrologic processes to nutrient cycling requires an understanding of 
biogeochemical transformations along hillslope flow pathways. Of primary importance is 
the impact of hillslope water residence time on oxidation - reduction (redox) conditions. 
Changes in the redox state of dissolved nutrients can result in losses of nutrients to the 


12 



atmosphere (e.g., DOC to CO 2 or NO 3 to N 2 0) or as storage (e.g., DOC to soil organic 
carbon or DOC to sorbed DOC). For example, Webster et al. (2008a, b) found that 
topographic controls on soil moisture and the quantity and quality of carbon substrates 
were important determinants of C0 2 efflux from soils. Transiently wet areas adjacent to 
wetlands, streams and lakes were a major source of soil CO 2 because of synchronicity in 
optimal temperature and moisture conditions during the growing season and a large pool 
of high quality substrate within freshly fallen leaves that accumulated in these 
depositional areas. In contrast, upland areas had insufficient soil moisture, which limits 
C0 2 production due to lack of soluble substrate; however, excess water in the wetland 
areas also limits C0 2 production due to lack of oxygen. Similarly, topography regulates 
N 2 0 efflux. During winter (Fairweather, 2007) and snowfree (Casson, 2008) periods of 
the year, the transiently wet areas adjacent to wetlands and the wetlands themselves were 
found to be the major source of soil N 2 0 efflux, and wetlands with redox conditions 
promoting denitrification processes were a major source of N 2 . 

Defining general rules that link hydrologic processes to nutrient export requires 
tools that precisely and accurately detect topographic characteristics within catchments. 
Digital terrain analysis generates terrain attribures related to hydrological processes that 
can assist in relating topographic composition and configuration to hydrological and 
biogeochemical response characteristics (Wagener et al., 2007). The potential of digital 
terrain analysis has been revolutionized by light detection and ranging (LiDAR), an 
airborne remote sensing technology that measures the distance to surfaces by illuminating 
the surface with laser beams and then analyzing the backscattered light. LiDAR has been 
used to produce fine-resolution digital elevation models of surfaces from which subtle 
topographic features can be detected. From LiDAR digital elevation models, terrain 
analysis can better define catchment boundaries, represent low order streams, delineate 
local depressions that form potentially wet areas, and identify ephemeral and permanent 
streams that connect these features to the stream (c.f., Creed and Sass, 201 1). For 
example, small hydrological features ( 1 st order streams, small wetlands) that were 
previously not detected, particularly in areas where vegetation is dense or the ground is 
otherwise shielded from aerial view (Bishop et al., 2008), can now be mapped. Thus, 
digital terrain analysis can now be used to derive topographic indicators for tracking 
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hydrological and biogeochemical dynamics across topographically simple to complex 
landscapes. Such knowledge will assist in generalizing process controls so that they can 
be applied to ungauged catchments (Tetzlaff el al . , 2008). 

Creed and Beall (2009) described the conceptual basis of topographic controls on 
NO3 -N and DON export in undisturbed headwater catchments on a sugar maple forested 
landscape. Topography influences (a) hydrological flushing potential (i.e., effect of 
topographic slopes on potential for variable source area to expand and tap into previously 
untapped nutrient source areas in the forest floor and surface mineral layers) and (b) 
hydrological storage potential (i.e., effect of topographic flats/depressions in catchment 
which represent wetlands that may either result in increased runoff when the water table 
in the wetland is near or at the surface, or decreased runoff when the water table in the 
wetland is well below the surface preventing water contributing to the wetland from 
being exported to the stream). This study builds on Creed and Beall (2009) by exploring 
similarity in topographic controls on other nutrient species, DOC and TDP, and 
determining if topographically regulated shifts in hydrological flushing versus storage 
potentails can explain differences in NO3 -N, DOC, DON and TDP export among 
catchments. The potential significance of this study is that it presents simple topographic 
indicators that can be used to track nutrient sources, sinks, and their transport and export 
to surface waters from both gauged and ungauged catchments on forest landscapes. 

2.2 Study Area 

The Turkey Lakes Watershed (TLW) is a 10.5 km” experimental watershed on the 
Algoma Highlands, about 60 km north of Sault Ste. Marie, Ontario (Figure 2.1). This 
watershed is located in the northern portion of the Great Lakes-St. Lawrence Forest 
Region, the second largest forest region in Canada, which extends from southeastern 
Manitoba to the Gaspe Peninsula. Climate is continental with mean annual precipitation 
and temperature of 1200 mm and 5.0 °C, respectively (Creed and Beall, 2009), and is 
strongly influenced by the proximity of Lake Superior to the west of the watershed and 
local orographic effects at locations of high relief. There is typically snowpack from late- 
November, early December through to late March, early April. Peak stream discharge 
occurs during snowmelt and again in October to November during autumn storms. 
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Figure 2. 1 : Map of the Turkey Lakes Watershed depicting the thirteen experimental catchments 
(from Creed and Beall, 2009). 

Bedrock geology is primarily Precambrian silicate greenstone, except north of the 
Batchawana Lake area and near the main inflow to Little Turkey Lake, where there are 
small outcrops of more felsic igneous rock (Semkin and Jeffries, 1983). Surficial geology 
of upland areas are glacial till with a depth < 1 m and frequent bedrock exposure; lowland 
areas are glacial till of 1-2 m on average; and bedrock depressions are characterized with 
deep till reaching up to 71 m in depth (Elliot, 1985). The glacial till is double-layered, 
comprising sandy loam ablation till overlying a compacted lower silt loam basal till. The 

'j 

upper layer is thin (< lm) and permeable (10' cm/sec), whereas the lower layer is thicker 
with lower permeability (10~ 5 cm/sec) (Johnston and Craig, 1986; Craig and Johnston, 
1988). The glacial till has been weathered to present-day loam, sandy loam and silty loam 
soil textures. The soils on the tills are Orthic Ferro-Humic and Humo-Ferric podzols with 
dispersed pockets of Ferric Humisols found in bedrock-controlled depressions and 
adjacent to streams and lakes (Wickware and Cowell, 1985). Forest is comprised of 
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mature to overmature sugar maple (>90%). Stand density (904 stems ha" 1 ), dominant 
height (20.5 m), diameter at breast height (15.3 cm), and mean basal area (25.1 m 2 ha" 1 ) 
are relatively uniform across vegetation types. The sparse understories of upland stands 
are dominated (>95%) by saplings and seedlings of sugar maple with a depauperate herb 
flora and a variety of ferns. The wetland understories are composed of the seedlings and 
saplings of the overstory trees, various ferns, herbs and a mix of feather and sphagnum 
mosses (Wickware and Cowell, 1985). 

The landscape exhibits substantial topographic variability, with gentle to steep 
hillslopes (11.5° to 20.8° in the experimental catchments) of different curvatures and 
lengths draining into streams. These hillslopes often drain through topographic flats 
and/or depressions containing organic -rich mineral or organic soils that are often 
transiently or permanently saturated. The presence of these topographic flats and 
depressions influences hydrological flow partitioning and pathways (Lindsay et al., 2004; 
Lindsay and Creed, 2005) and biogeochemical cycling (Webster et al., 2008a, b). 

The TLW has been used as an experimental watershed by the federal government 
agencies since 1980, initially to study the potential impact of acid rain and climate 
change in terrestrial and aquatic ecosystems in the Canadian Shield area (Foster, 1985; 
Jeffries et al., 1988). Since 1980, daily meteorological data (maximum and minimum 
temperature, precipitation and solar radiation data) around the watershed have 
continuously been measured at the meteorological station close to the watershed. In 
addition, since 1981, discharge and various water chemistry parameters have been 
monitored for 13 gauged headwater catchments within the TLW (Figure 2.1). A 
harvesting experiment was carried out in 1997 at three of the headwater catchments (c31, 
c33 and c34) to investigate the hydrological consequences of different forest management 
alternatives: a diameter limit harvest, shelterwood harvest and a selection harvest. For 
this study, pre-harvest data from 1981-1997 are used for c31, c33 and c34, while data 
from 1981-2008 are used for the ten, unharvested catchments. 
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2.3 Methods 


2.3.1 Topographic Indicators 

A 5 m Light Detection and Ranging (LiDAR), hydrologically conditioned digital 
elevation model (DEM) with vertical accuracy of 0.15 m in open canopy and 0.30 m in 
closed canopy was used to derive topographic metrics. Pits and depressions were filled 
(Planchon and Darboux, 2002), and the D8 algorithm was used to determine catchment 
contributing areas and surface drainage networks (O’Callaghan and Mark, 1984; Jenson 
and Domingue, 1988). Specific contributing area of 2000 nr and stream length threshold 
of 25 grid cells (125 m) were selected because they resulted in drainage networks that 
matched most closely those observed in the catchments. 

Topographic indicators representing hydrological storage mechanisms (i.e., 
topographic flats/depressions that were assumed to be wetlands) were identified 
following a stochastic probability approach (Lindsay and Creed, 2005) (Figure 2.2). The 
approach involved deriving a grid map of the probability of depression (pdep ) using 
Monte Carlo simulation where (1) the elevation error term was stochastically drawn from 
a frequency distribution of elevation error terms with a standard deviation equal to the 
vertical accuracy of the DEM (i.e., 0.30 m); (2) this error term was added to the DEM; 

(3) the flats/depressions in this modified DEM were filled and each grid cell modified by 
the depression filling process identified; (4) a probability of depression occurrence (pdep ) 
was calculated by the number of times each grid cell was contained in a depression, with 
pdeps ranging from 0 (area with no probability of being a depression) to 1 (area with 
100% 
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Stream Wetland 


Figure 2.2: Method used to derive topographic metrics for hydrologic storage potential. 
Random errors were introduced into the DEM to calculate depressions, this process was 
repeated iteratively and the probability of depression was calculated, and wetlands were 
delineated using a 0.3 threshold determined by the minimum difference in area between 
simulated and observed (ground truthed) wetlands. 
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probability of being a depression). The Monte Carlo simulation was repeated using 
different elevation error terms until a stable solution occurred, defined by when the RMS 
difference in pdep between two consecutive simulations was < 0.001. Grid cells with 
pdep values greater than a defined critical threshold were identified as topographic 
depression and/or flat areas (e.g., wetlands). The critical threshold in pdep (0.30), above 
which a depression was defined, was determined by calculating the lowest area difference 
and highest Kappa statistics between stochastically derived and field surveyed wetlands. 

Topographic indicators representing hydrological flushing mechanisms were 
derived using the following methods (Figure 2.3). Topographic Index (TI) is an indicator 
of topographically driven soil moisture conditions and was calculated as ln(a/tan//), 
where a is upslope contributing area of a given site based on drainage directions 
calculated using the D-infinity algorithm (Tarboton, 1997) from a specific contributing 
area and slope gradient maps and (i is local slope angle (Beven and Kirkby, 1979). The TI 
grid map of each catchment was normalized (TI n ) by subtracting catchment average TI 
values from specific TI for each grid cell. Variable Source Area (VSA) is an indicator of 
potential hydrological flushing areas were estimated by starting at the stream grid cells, 
where the TI n values are the highest, and recursively moving outward from the stream 
grids to grid cells of lower TI n value until a point where the TI n values start to increase 
again or the catchment boundary is reached (Figure 2.4A). The proportion of VSA where 
hydrological flushing is likely to occur on an annual basis was referred to as effective 
variable source area (effVSA) and was estimated by grid cells with TI n values > the 75 th 
percentile of the TI n frequency distribution within the VSA, a percentile that was 
previously found to provide the best prediction for NO 3 -N export (Creed and Beall, 

2009) (Figure 2.4B). A situation was possible where the same pixel can function as both 
hydrological storage and hydrological flushing (i.e., if topographic flats/depressions 
occurred within the effVSA). For this reason, both effVSA with and effVSA with no 
topographic flats/depressions (i.e., no wetlands) (NW<?//VSA) were derived and 
considered in the development of models (Figure 2.4C). The potential rate of effVSA 
expansion and contraction (d 2 <?//VSA/dTI n 2 ) was calculated as a second derivative of the 
best-fit polynomial to the TI n frequency distribution within the effVSA. 
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Figure 2.3: Description of methods used to derive topographic metrics for hydrological 
flushing potential: (A) VSA is calculated spatially from the TI n map by starting at the 
stream and going recursively to lower Tin values until there is a breakpoint where values 
increase; (B) effVSA is the upper quartile of the TI n values within the VSA and describes 
the portion of the VSA where return flow is likely to occur on an annual basis; 
NW<?//VSA is calculated spatially as the area of the effVSA outside wetlands; and (C) the 
rate of change of effVSA is calculated as the second derivative of the polynomial best fit 
of curve. 
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Figure 2.4: Spatial variation in annual Nitrate-N, dissolved organic carbon (DOC) 
dissolved organic nitrogen (DON), and total dissolved phosphorus (TDP) in the Turkey 
Lakes Watershed. 
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Ratios of the topographic metrics were also considered because the relative 
presence of one metric with respect to another might increase predictive capacity. For 
example, the ratio of topographic sinks of a nutrient (e.g., proportion of topographic 
flats/depressions that denitrify NOf-N) to topographic sources of a nutrient (e.g., 
proportion of VS A, effVSA or NWc//VSA that are sources of NO3 -N) were derived to 
see if they improved the predictive capacity of topographic metrics for NOa'-N export. 

All digital terrain analyses were performed using Terrain Analysis System (TAS) 
software (Lindsay, 2005). 

2.3.2 Solute Export 

The concentration of nutrients discharged into the streams were determined from samples 
collected at mid-day every two weeks during the winter, daily during spring snowmelt, 
and weekly or every two weeks during the summer and autumn. Samples were collected 
at the same sampling point, the center of the stream, for each stream. From each sample, 
particulate matter was removed by filtration through Fisher Q 8 (coarse, fast flow) paper 
filters. Samples were analyzed within 48 h of collection for NH 4 + -N, NOf-N. and total 
dissolved nitrogen (TDN) using sodium nitroprusside, cadmium reduction and cadmium 
reduction after autoclave digestion methods, respectively, on Technicon Autoanalysers. 
DON was calculated from TDN minus dissolved inorganic nitrogen (DIN = NFL + -N + 
NO 3 -N). DOC was determined by removing dissolved inorganic carbon (DIC) by 
purging with N 2 after acidification, converting DOC to DIC by persulfate oxidation 
catalyzed by UV, then converting the resulting DIC to CO 2 by acidification, which was 
measured by colorimetry. TDP was analysed on a Technicon Autoanalyser IIC after 
autoclave digestion using the molybdophosphoric acid colour reaction. 

Daily concentrations of nutrients were estimated by interpolating measured values 
directly before and after the unknown value. Daily fluxes of nutrients were then estimated 
by the product of the total daily discharge from continuously logged stream stage at V- 
notch weirs on the catchments and the interpolated daily concentrations of nutrients. 
Annual nutrient fluxes of nutrients were calculated as the sum of daily flux for the water 
year (June 1 to May 31) from 1981 to 2008 (except for the three harvested catchments, 
c31, c33 and c34, where the period was from 1981 to 1996 water years). 
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2.3.3 Relationships between Topographic Indicators and Nutrient Export 

The nature, strength and significance of the relationships between the catchment’s 
topographic indicators and its nutrient export characteristics were evaluated using 
correlations and linear regressions (p = 0.05) using Sigma Plot (Version 11, Systat 
Software Inc., 2008). 

2.4 Results 

2.4.1 Topographic Indicators 

Topographic indicators of hydrological storage potential (i.e., topographic flats/ 
depressions) varied substantially (Table 2.1). Catchments c37, c38 and c50 had 
substantially more wetlands (>10% topographic flats/depressions), whereas c33 and c47 
had very few wetlands (<1% topographic flats/depressions). The distribution of wetlands 
within the catchments varied, with some containing large contiguous wetlands (e.g., c38) 
and others smaller discontiguous wetlands that cascaded one to the next via the surface 
drainage network (e.g., c50). 

Topographic indicators of hydrological flushing varied substantially (Table 2.1). 
VSA was smallest in c47 (23.62%) and largest in c31 (61.66%). The patterns in ejfV SA 
did not follow VSA, reflecting differences in concavity versus convexity of the lower 
hillslopes draining into the stream. The effVSA was smallest in c38 (6.34) and c37 
(6.82%), where water moves down steep hillslopes, draining into much more gentle 
hillslopes before exiting the catchment, and largest in c33 (16.04%), where steep to 
moderate slopes drain directly to the stream. Topographic flats/depressions that would 
detain the movement of water from the hillslopes to the stream were prevalent in some of 
the effVSAs. When these topographic flats/depressions were removed from effVSA, the 
proportion of area contributing variable source area runoff decreased among all 
catchments (except c47 where it stayed the same), but the rank order of the NW<?//VSA in 
catchments generally remained the same, with c37 (1.73%) and c38 (1.55%) having 
lowest NWfjOWSA and c33 having highest NW<?//VSA (15.75%). These topographic 
indicators of VSA were correlated with one another (Table 2.2), because they were 
refinements of the same process. For this reason they were not combined in the models 
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Table 2.1: Topographic metrics and ratios or topographic metrics among headwater 
catchments in the Turkey Lakes Watershed (VS A = variable source area; effV SA = 
effective VS A; NW<?//VSA = area of effVSA without wetland; d 2 <?//VSA/dTI 2 n = rate of 
effVSA expansion and contraction). 


Catchment 

Wetlands 

(%) 

VSA 

(%) 

effVSA 

(%) 

NWeffVSA 
(%) ' 

W: 

VSA 

W: 

NeffVSA 

dffffVSA/dTlff 

c31 

2.88 

61.66 

13.10 

11.33 

0.05 

0.22 

2.00E-06 

c32 

1.00 

54.54 

14.32 

14.05 

0.02 

0.07 

4.00E-06 

c33 

0.50 

45.13 

16.04 

15.75 

0.01 

0.03 

2.00E-05 

c34 

1.12 

36.04 

9.70 

9.22 

0.03 

0.12 

4.00E-05 

c35 

1.06 

36.38 

13.18 

13.12 

0.03 

0.08 

4.00E-06 

c37 

15.01 

33.63 

6.82 

1.73 

0.45 

2.20 

1.20E-06 

c38 

20.54 

31.57 

6.34 

1.55 

0.65 

3.24 

6.00E-07 

c39 

5.97 

47.13 

12.96 

11.05 

0.13 

0.46 

4.00E-06 

c42 

8.48 

51.84 

13.71 

10.06 

0.16 

0.62 

4.00E-06 

c46 

1.35 

30.57 

9.73 

9.08 

0.04 

0.14 

8.00E-06 

c47 

0.36 

23.62 

7.58 

7.58 

0.02 

0.05 

2.00E-06 

c49 

3.97 

47.27 

9.65 

8.62 

0.08 

0.41 

2.00E-06 

c50 

10.03 

38.09 

7.84 

3.72 

0.26 

1.28 

1.00E-06 

Average 

5.56 

41.34 

10.84 

8.99 

0.15 

0.69 

3.08E-06 

Max 

20.54 

61.66 

16.04 

15.75 

0.01 

0.05 

4.00E-05 

Min 

0.36 

23.62 

7.58 

1.55 

0.65 

3.24 

6.00E-07 
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Table 2.2: Pearson correlation between average topographic metrics, including ratios between metrics, for 13 catchments in the 
Turkey Lakes Watershed (VSA = variable source area; effVSA = effective VSA; NW effVSA = area of effV SA without wetland; 
d 2 <?//VSA/dTI 2 n = rate of effVSA expansion and contraction; NS = not significant, all other relationships p < 0.05). 


Topographic Metrics 

VSA (%) 

effVSA (%) 

NWc/fVSA (%) 

W: VSA 

W: effVSA 

W: NWcifVSA 

A 1 effV S A/dTI n " 

Wetlands (%) 

NS 

-0.575 

-0.813 

0.986 

0.978 

0.932 

NS 

VSA (%) 


0.707 

NS 

NS 

NS 

NS 

NS 

effVSA (%) 



0.929 

-0.628 

-0.644 

-0.629 

NS 

NWefjfVSA (%) 




-0.828 

-0.834 

-0.787 

NS 

W:VSA 





0.998 

0.977 

NS 

W:effVSA 






0.983 

NS 

W:NWe#VSA 







NS 
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constructed to predict nutrient export. Topographic properties of hillslopes that would 
affect the rate of potential expansion of effV SA (d 2 e//VSA/dTI n 2 ) varied by orders of 
magnitude, from high rates of potential expansion in c33 (2 x 10" 5 ) and c34 (4 x 10 5 ) 
(tapping into relatively high NO 3 -N generating areas that exist on less concave (less flat) 
hillslopes) to low rates in c38 (6 x 10" ) (tapping into relatively low NCf-f-N generating 
areas that exist on more concave (more flat) hillslopes) (Figure 2.2). 

A broad range existed in the topographic indicators of hydrological flushing 
versus storage potential among the catchments. Catchments c38, c37 and c50 had the 
largest proportion of wetlands, but the smallest effVSA and HW effVSA, and the least 
expandable effVSA (d 2 c//VSA/dTI n 2 ). For example, c38 had the largest wetland (20.54%) 
but among the smallest VSAs (31.6%), the smallest effVSA (6.3 %) and the smallest 
NWeffVSA (1.55%). In contrast, c33 had among the smallest wetlands (0.5%) but the 
largest VSA (45.1%), the largest effVSA (16.0%), giving it the largest NWejfVSA 
(15.8%), and one of the most expandable effV SA (d 2 c//VSA/dTI n 2 ). Catchments c33 and 
c38 represent extremes in terms of topographic indicators of hydrologic flow, and 
therefore comparisons of these two sites are highlighted in subsequent analyses. 

When the ratio of the proportion of topographic flats/depressions (where rising 
water table leads to inundation) to proportion of VSAs (where rising water table leads to 
expansion into unsaturated soils) (W:effVSA) was calculated the magnitude of difference 
among the catchments increased. For example the W: effVSA for c33 was 0.03 
(emphasizing the role of effV SA relative to W); whereas, the W:effVSA for c38 was 3.24 
(emphasizing role of W to effVSA). 

2.4.2 Nutrient Export 

Nutrient export from the catchments varied substantially (Table 2.3). Catchments 
representing lows in inorganic nutrient export tended to be the highest in organic nutrient 
export, and vice versa. For example, catchment NO 3 -N export ranged from a low of 1.02 
kg ha " 1 a " 1 (c38) to a high of 5.15 kg ha " 1 a " 1 (c35). In contrast, catchment organic 
nutrient export ranged from a low of 13.19 kg ha " 1 a " 1 for DOC, 1.02 kg ha " 1 a " 1 for DON, 
and 0.0169 kg ha 1 a 1 for TDP in catchment c35, to a high of 47.22 kg ha 1 a 1 for DOC, 
2.00 kg ha " 1 a " 1 for DON, and 0.0558 kg ha 1 a 1 for TDP for c38 (Table 2.3). 
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Table 2.3: Annual water and nutrient export (coefficient of variation) of dissolved 
organic carbon (DOC), dissolved nitrate (NO 3 -N), dissolved organic nitrogen (DON) and 
total dissolved phosphorus (TDP) from the Turkey Lakes Watershed for water years 
(June through May) from 1981 to 2008, except for catchment c31, c33 and c34 where 
data from 1981 to 1996 were used due to disturbance. 


Catchment 

Area 

Water 

NCL-N 

DOC 

DON 

TDP 


(ha) 

(mm a' 1 ) 

(kg ha 1 a 1 ) 

(kg ha 1 a 1 ) 

(kg ha 1 a 1 ) 

(kg ha 1 a' 1 ) 

c31 

4.94 

569 (21) 

3.98 (26) 

13.51 (24) 

1.06 (29) 

0.0174 (37) 

c32 

6.5 

495 (30) 

4.56 (35) 

10.00 (25) 

0.77 (55) 

0.0138 (49) 

c33 

23.38 

531 (19) 

5.19(20) 

13.13(15) 

1.03 (28) 

0.0143 (48) 

c34 

68.59 

663 (14) 

5.43 (19) 

14.16(14) 

1.26 (49) 

0.0186 (38) 

c35 

4.02 

632 (26) 

5.15(29) 

13.19(29) 

1.02 (41) 

0.0169 (54) 

c37 

15.36 

633 (21) 

1.88 (33) 

34.11 (23) 

1.60 (26) 

0.0340 (25) 

c38 

6.46 

612 (26) 

1.02 (43) 

47.22 (24) 

2.00 (29) 

0.0558 (30) 

c39 

17.25 

518 (24) 

3.93 (30) 

15.92 (24) 

0.87 (35) 

0.0165 (31) 

c42 

18.52 

457 (22) 

2.84 (41) 

20.61 (22) 

0.96 (33) 

0.0165 (26) 

c46 

43.19 

654 (29) 

3.18(36) 

14.22 (41) 

1.09 (49) 

0.0199 (40) 

c47 

3.43 

555 (26) 

2.72 (38) 

12.72 (32) 

0.83 (37) 

0.0135 (40) 

c49 

14.81 

710(29) 

2.93 (39) 

18.90 (42) 

1.16(34) 

0.0228 (31) 

c50 

9.47 

776 (22) 

2.02 (40) 

32.48 (19) 

1.67 (31) 

0.0344 (30) 

Average 

18.1 

600 

3.45 

20.01 

1.18 

0.0227 

Maximum 

68.59 

776 

5.43 

47.22 

2.00 

0.0558 

Minimum 

3.43 

457 

1.02 

10.01 

0.77 

0.0135 
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Catchment differences in nutrient export were also depicted taking the residuals 
between each catchment and the average across all catchments for each nutrient (Figure 
2.4). For example, c35 had 49% higher NO f-N export than the average export among 
catchments, while c38 was 70% lower in NO3 -N export than the average. The higher 
elevation catchments (e.g., c42, c46, c47, c49) that had little to no topographic 
flats/depressions were the exception, showing little difference from average export or 
slightly lower export than the average export across all catchments (Figure 2.5). NOf-N 
export behaviour was distinct from DOC, DON and TDP export behaviour. For example, 
c35, which had higher NO.F-N export, had lower export of organic nutrients. Spatially, 
c35 exported 14% less DON, 34% less DOC, and 26% less TDP than the average export 
across all catchments; whereas, c38 exported 69% more DON, 136% more DOC and 
146% more TDP than the average export across all catchments (Figure 2.5). 

2.4.3 Topographic Indicators vs. Nutrient Export 

There were generally no significant correlations between topographic indicators and 
water export (Table 2.4); the one exception was effNSA versus water export, where 
increase in effNSA led to a decrease in water export (r 2 =0.379, p<0.05), but an increase 
in HW effNSA had no significant effect on water export. However, there were significant 
correlations between topographic indicators and nutrient export (Table 2.5). For example, 
for NO.f-N export, the best linear model was the simple indicator of NW effNSA (r = 
0.759, p < 0.001), while the best multiple linear model was the combination of simple 
indicators of NWcjffVSA and d 2 e//VSA/dTI n 2 (R 2 = 0.907, p < 0.001). Compound 
indicators, including W:VSA or N/: effNSA. did not improve the linear or multiple linear 
models (Table 2.5). For DOC, DON and TDP export, the best linear model using simple 
indicators was topographic flats/depressions within the catchment (R = 0.935 for DOC, 
0.694 for DON, and 0.821 for TDP). Multiple linear models did not improve performance 
However, compound indicators W:VSA or W:c//VSA did improve performance. For 
DOC, DON and TDP export, the best linear model using compound indicators was 
W:effNSA (R 2 = 0.957 for DOC, 0.783 for DON, and 0.912 for TDP), which resulted in 
improvements of up to 10% in the explanation of variation in organic nutrient export 
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Average Nutrient Export (kg/ha/yr) 



o DOC = -2.2x + 39.3, p<0.001, r 2 = 0.746 

o N0 3 -N = 0.3x + 1 .0, p<0.001 , r 2 = 0.758 

• DON = -O.lx + 1.8, p<0.001, f = 0.684 

o TDP = - 0.002x + 0.043, p<0.001 , r 2 = 0.688 

• Water 


Figure 2.5: Average water yield and nutrient export versus No Wetland effective Variable 
Source Area (NWejffVSA) in 13 catchments of the Turkey Lakes Watershed. 
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Table 2.4: Relationships between topographic metrics and annualized water export from 
headwater catchments in the Turkey Lakes Watershed (Direction = direction of 
relationship, r = coefficient of determination, and NS = not significant, all other 
relationships p < 0.05). 


Topographic Metrics 

Water 


Direction 

2 

r 

Wetlands (%) 

NS 

NS 

VS A (%) 

NS 

NS 

effVSA (%) 

- 

0.379 

NWeffVSA (%) 

NS 

NS 

d 2 <?f/VSA/dTI n 2 

NS 

NS 

W:VSA 

NS 

NS 

W:effWSA 

NS 

NS 

NWeffVSA and dT//VSA/dTI n 2 

NS 

NS 
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Table 2.5: Relationships between topographic metrics and annualized nutrient export 
from headwater catchments in the Turkey Lakes Watershed (Dir. = direction of 
relationship, r = coefficient of determination, and NS = not significant, all other 
relationships p < 0.05). 


Topographic 

Metrics 

no 3 '-n 

DOC 

DON 

TDP 


Dir. 

r 

Dir. 

r 

Dir. 


Dir. 

r 

Wetlands (%) 

- 

0.646 

+ 

0.935 

+ 

0.694 

+ 

0.821 

VSA (%) 

NS 

NS 

NS 

NS 

NS 

NS 

NS 

NS 

effVSA (%) 

+ 

0.561 

- 

0.456 

- 

0.489 

- 

0.484 

NWeffV S A (%) 

+ 

0.759 

- 

0.747 

- 

0.682 

- 

0.687 

d 2 c#VSA/dTI n 2 

+ 

0.394 

NS 

NS 

NS 

NS 

NS 

NS 

W:VSA 

- 

0.619 

+ 

0.951 

+ 

0.752 

+ 

0.887 

W'.ejfVSA 

- 

0.618 

+ 

0.957 

+ 

0.783 

+ 

0.912 

NWeffV S A and 
dreffV S A/dTI n 2 

+ 

0.907 

NS 

NS 

NS 

NS 

NS 

NS 
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(Table 2.5). Among the organic nutrient exports, model results for DON explained less 
variation (about 13-17% less) than DOC and TDP (Table 2.6). 

2.5 Discussion 

Catchment-based water quality monitoring programs are few and far between, and 
increasingly at risk of being closed due to budget constraints. While perceived as 
“expensive and wasteful” (Lovett et al., 2007) these long-term monitoring programs are 
in fact unique and valuable scientific assets, providing insights into spatio-temporal 
hydrological and biogeochemical dynamics that cannot be gained using other research 
approaches, such as experimentation or modeling. Importantly, the value and impact of 
these scientific assets can be increased by developing a “similarity metric” approach in 
that heterogeneities in environmental factors (e.g., topography, geology and climate) are 
related to catchment functional responses (Tetzlaff et al., 2008). By developing a metric- 
based approach we not only capitalize on the unique long-term records these studies 
provide, but also provide a means to extract general principles that can be used in 
ungauged basins, thereby increasing efficiencies by removing the need to expand or 
create new monitoring programs. 

Using catchment classification and similarity for predictions in ungauged 
catchments based on the characterization of appropriate metrics of both temporally 
variable (e.g., timing of extreme precipitation events, Castellarin et al., 2001) and 
spatially heterogeneous (e.g., soil and topography, Larsen et al., 1994) catchment 
properties on runoff responses have a long tradition in hydrology. The strength of this 
approach has improved with advances in DEM resolution and new digital terrain analysis 
methods (e.g., Lindsay and Creed, 2005; Creed and Beall, 2009) allowing for precise 
characterization of the topography used in these metrics. Although considerable effort has 
been made to apply metrics for defining similarities in hydrologic response (Larsen et al., 
1994; Castellarin et al., 2001), little work has been done on developing metrics relevant 
for predicting nutrient export dynamics, despite various studies relating nutrient export to 
environmental controls of spatial and temporal nature (Heathwaite and Johnes, 1995; 
Fraterrigo and Downing, 2008; Andrews et al., 2011). 
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Table 2.6: The best regression models for predicting average annual nutrient export from 
headwater catchments in the Turkey Lakes Watershed. All models were significant at 
pcO.OOl. 


Response Variable 

r 

Equation 

NCL-N (kg N ha 1 a 1 ) 

0.907 

NO 3 --N = 0.978 + (0.235 x NWeffW SA) + 
(50557.598 x d 2 c//VSA/dTl 2 n ) 

DOC (kg C ha ' 1 a 1 ) 

0.957 

DOC = 12.499 + (10.956 x W: effV SA) 

DON (kg N ha 1 a' 1 ) 

0.783 

DON = 0.954 + (0.326 x W: effVSA) 

TDP (kg P ha 1 a' 1 ) 

0.911 

TDP = 0.0146 + (0.01 17 x W: effVSA) 
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We examined the effect that topography has on nutrient export from headwater 
catchments in order to find key metrics that can be applied to prediction in ungauged 
basins. We found considerable variation in nutrient export at the TLW, despite 
uniformity in climate, geology, forest type and forest age (Creed et al., 2002). 
Topographic characteristics were largely responsible for explaining variation in nutrient 
export among catchments, with topography controlling hydrological processes that 
influence the production, availability and delivery of nutrients to the downstream surface 
waters. In particular metrics related to hydrological storage potential (i.e., topographic 
flats/depressions) and hydrological flushing potential (i.e., NW<?//VSA and 
d“<?//VSA/dTI n ) explained the majority of variation in nutrient export. 

Wetlands that form in topographic flats/depressions have been identified as a sink 
of NC> 3 ~-N and a source of DOC, DON and TDP. Previous studies have shown that 
wetlands can act as sources (Lee et al., 1975; Sloey et al., 1978), sinks (Peterjohn and 
Cornell, 1984; Kadlec, 1986) and/or transformers of nutrients. Their importance in 
nutrient export is disproportionate to their area on the landscape. Although wetlands 
commonly only cover a small proportion of the landscape they are usually positioned so 
that much of the discharge waters and the associated nutrients pass through them on their 
way to the receiving water bodies (Devito et al., 1989). The wetland in c38 from our 
study area was identified as a sink of NCff-N. suggesting that NOf-N from upland 
drainage waters that is intercepted in the wetland is transformed and subsequently 
released via denitrification to NoO and/or N 2 (Creed, Unpubl. Data). Alternatively, NO 3 - 
N from upland drainage waters that is intercepted in the wetland is retained via abiotic 
immobilization (Davidson et al., 2003) and/or biotic immobili z ation (Gundersen and 
Bashkin, 1994; Sirivedhin and Gray, 2005), although we did not quantify the importance 
of these mechanisms of nitrogen retention. 

In contrast, wetlands were a source of DOC, DON and TDP, with the magnitude 
of nutrient export proportional to wetland area. There are several mechanims to explain 
this observation. The physical and chemical environment within wetlands affects the 
solubility, lability and mobility of nutrients. Of particular importance are the organic 
deposits as a major source of soluble organic by-products (Qualls and Richardson, 2003). 
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For example, saturation may suppress mineralization of DON to inorganic nitrogen 
(NH 3 + -N) (Yang and Fan, 2003) leading to organic forms accumulating in wetlands. The 
wet conditions also have a large impact on the redox state of nutrients (Sorrell and 
Armstrong, 1994; Reddy and DeLaune, 2008). Low redox potentials may also facilitate 
abiotic conversion of N0 3 "-N to organic forms of nitrogen through the ferrous wheel 
hypothesis, which involves reduction of Fe (III) hydroxide by organic matter in forest 
floor to release Fe (II) (Davidson et al., 2003). Low redox potentials may also facilitiate 
biotic conversion of N0 3 "-N to gaseous forms of N 3 0 and N 3 . The potential for 
alternative fates of nitrogen from uplands passing through wetlands en route to the stream 
may explain the relatively weak relationship between wetland proportion and DON 
compared to DOC and TDP. 

Uplands and their VS As (Hewlett and Hibbert, 1967) have been important in 
predicting nutrient export in many studies (e.g., Govindaraj, 1996; Walter et al., 2000; 
Marjerison et al., 201 1). In this study, despite noticeable variation in VSA proportion 
among the studied catchments, VSA proportion did not explain the variation of nutrient 
export among catchments. However, we did find that different derivatives of VSA were 
important indicators in predicting nutrient export. 

For example, the proportion of effV SA, which represents the proportion of 
frequently flushed area within the VSA, was more effective in predicting N0 3 "-N export. 
The direct relationship observed between effVSA and NOf-N export suggested an 
increased rate of hydrological flushing of N0 3 '-N from catchments with larger effVSAs. 
Aerobic conditions in effVSA surfaces, when not saturated, favor nitrification (oxidation 
of NH 3 + -N to N0 3 ~-N), making effVSA a net source of N0 3 "-N. Nitrate accumulates in 
these areas until sufficient precipitation events cause the water table to rise and flush the 
accumulated N0 3 ’-N for export to streams (Creed et al. 1996; Creed and Band, 1998). If 
not absorbed by plants, N0 3 "-N, which naturally possesses less affinity for adsorption to 
mineral soils (Foster et al., 1982; Andrews et al., 2007), can be easily mobilized from 
effVSA to receiving water bodies. 

Furthermore, the proportion of effVSA without topographic flats/depressions 
(NW effVSA) was even more effective in predicting NO f-N export. EffVSA may contain 
wetlands, and NWeffVSA removed the potential effects of wetlands within effVSA. 
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HW effVSA explained more variation than either of these indicators did on their own. 
Coupling NW<?//VSA with potential lateral expansion of ef'fWSA (d 2 <?//' VSA/dTI n 2 ) 
further improved predictions of NO 3 -N export. The size of hydrologically connected 
portion of VSA ( effVSA ) varies with season and storm characteristics. Expansion of 
effVSA in response to climatic forcing can result in larger flushable areas and longer 
flushing time, causing greater NO 3 -N export (Creed and Band, 1998b). The role of 
expanding saturated VSA in terms of increasing flushable nutrient concentrations has 
also been reported for NO 3 -N by Sickman (2001) and Oeung el al. (2010) and for DOC 
by Boyer et al. (1995). 

It is important to recognize that the VSA and its derivatives are not related to 
water export (i.e., expansion of rising water table into upland areas does not correspond 
to an increase in water at the stream). Rather, as VSA increases, the amount of water 
exported from the catchment decreases as the water enters uplands areas that were 
previously dry and flushes NO 3 ' -N to the stream. Consequently, as VSA increases, water 
export tends to decrease while NOa'-N export increases. 

The different derivatives of VSA were also significantly related to organic 
nutrient and TDP export, but negatively correlated, suggesting different controls on their 
export. Generally catchments with large effVSA proportion contain no to few wetlands 
whereas those with small effVSA proportion contain substantial wetlands. The inverse 
relationship between effVSA and DOC, DON and TDP relative to N 03 ~-N may be linked 
to physical processes. For example, the increased adsorption of dissolved organic matter 
by mineral soils in upland dominated catchments may decrease exports of organic forms 
of nutrients from these catchments (McDowell and Likens, 1988; Fell man et al., 2009). 
As effVSA proportion increases, the potential of dissolved organic matter in upland 
dominated catchments to reach stream networks would decrease because the transport 
pathways are through mineral-rich soils that strongly adsorb dissolved organic matter 
(Wallis, 1979, Qualls et al., 2000). Mineral soil adsorption has been well recognized as a 
primary process responsible for reducing dissolved organic matter export (McDowell and 
Likens, 1988; Moore et al., 1992; Kalbitz et al., 2000; Moore and Turunen, 2004), while 
organic soils with little DOC adsorption capacity and the potential to generate soluble 
carbon compounds are known to be large exporters of DOM (Aitkenhead and McDowell 
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2000; Aitkenhead-Peterson et al., 2007). For example, due to efficient adsorption by 
mineral soils, percolation of DOC through subsurface soil can cause a typical decline of 
DOC concentration by 50-90% (McDowell and Wood, 1984; Neff and Asner, 2001). 

The inverse relationship between effV SA and DOC, DON and TDP relative to NOf-N 
may also be linked to biological processes. For example, aerobic conditions in effV SA 
soils of uplands would faciliate microbial organic matter mineralization providing an 
additional source of NO3 -N for export while depleting dissolved organic matter reserves 
faster than what is generated by decomposing detritus (Aiken et al., 201 1). 

However, wetland proportion was a more important metric than VS A derivatives 
for predicting DOC, DON and TDP export. We developed new metrics reflecting the 
relative proportion of wetland area with respect to VSA (W:VSA) and effV SA 
(W:effVSA). In general, W:VSA and W:<?//VSA metrics were more effective in 
predicting export of organic forms of nutrients (DOC, DON, TDP) than wetland 
proportions (% wetland). While the improvement was marginal in DOC prediction (-2%), 
it was larger for DON and TDP with W:VSA and W:effVSA, improving predictions by 
~6% and 9%, respectively. Catchments with larger wetland area with respect to their 
VSA and effV S A tend to export more organic nutrients; the sources have a greater 
influence than the sinks of DOC, DON and TDP. 

2.6 Conclusions 

We need tools to scale our process-based understanding from experimental catchments to 
the landscape where policy and managaement decisions are made. Topography is an 
important determinant of nutrient export through its regulation of hydrological processes 
that in turn influence formation of nutrient hot spots (sources for export) versus cold 
spots (sinks for export), transformation of the nutrients to dissolved and/or gaseous forms, 
and transport to the stream. Specifically, topography exerts an influence on hydrological 
storage potential (i.e., topographic flats/depressions where wetlands form allowing the 
accumulation of nutrient-rich organic deposits), and hydrological flushing potential (i.e., 
topographic controls on the effVSA and its potential for expansion where rising water 
table allows the prolonged flushing of nutrient-rich near surface/surface soils as rising 
water table reaches new soils). Topographic characteristics of individual catchments 
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create different potentials for hydrological storage versus flushing which contributes to 
the substantial variation in nutrient export among catchments, and the relative proportion 
of one with respect to the other should be considered to improve prediction of nutrient 
export. Recent developments have revolutionized topographic models of catchments 
which means that tools can be developed that incorporate topographic indicators 
representing hydrological and biogeochemical controls on nutrient export behavior of 
catchments that can be used as a basis for characterizing responses in ungauged 
catchments. 
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Chapter 3: Russian nesting dolls effect - Using wavelet 
analysis to reveal non-stationary and nested stationary 
signals in water yield from catchments in a northern 
forested landscape 1 


The previous chapter focused on the role of topography on spatial variation in long-term 
average annual nutrient export, whereas this chapter focuses on the role of topography in 
influencing temporal variation in annual water yield (as a precursor for examining 
annual nutrient export). This chapter tests the hypothesis that catchment water yields 
contain a complex mixture of non-stationary and stationary signals driven by climate that 
are manifested differently in catchments with different topographic characteristics. A 
novel analytical framework is introduced to deconstruct time series. Annual time series 
from catchments with topographic characteristics that represent the range in 
hydrological flushing potential, hydrological storage potential, and hydrological loading 
potential are then deconstructed into non-stationary (trend) and stationary (cycle) 
signals, and then climate drivers are related to these deconstructed signals to determine 
their contribution to these temporal variation in water yields. The significance of this 
chapter is that it presents a novel analytical framework for deconstructing export signals, 
so that potential for differential sensitivity of catchment water yields to climatic forcing 
can be established. 

3.1 Introduction 

Understanding climate effects on water resources continue to be a major agenda among 
hydrologists, ecologists and policy makers. Catchments in various regions have been 
affected by non-stationary climate warming linked to anthropogenic activities (Amell, 
2004; Christensen et al., 2004; Chen et al., 2006; Pike et al., 2008; Zwiers et al., 201 1). 
Catchments have also been affected by stationary climate patterns, driven by large-scale 


*A version of this chapter has been published (Mengistu SG, Creed IF, Kulperger RJ, Quick CG. 2012. 
Russian nesting dolls effect - Using wavelet analysis to reveal non-stationary and nested stationary signals 
in water yield from catchments on a northern forested landscape. Hydrological Processes ) 


45 




climatic oscillations caused by fluctuations in sea surface temperature or sea level 
pressure (Keener et al., 2010; Kahya, 2011; Niedzielski, 201 1) [e.g., Multivariate El Nino 
Southern Oscillation (ENSO) Index (MEI), Atlantic Multidecadal Oscillation (AMO), 
Northern Atlantic Oscillation (NAO), and Pacific Decadal Oscillation (PDO)], although 
there is some evidence that these oscillations may not be strictly stationary [e.g., ENSO 
(Gaucherel, 2010)]. Observed temporal variability in water yields may reflect the 
combined effect of non-stationary climate trends and stationary climate oscillations with 
multiple periodicities (Zhou et al., 2008). There is a critical need to discriminate between 
non-stationary and stationary signals in changing water yields from these catchments. 

Discriminating between non-stationary responses (deterministic responses where 
the statistical mean and variance change with time, predictably and unpredictably) and 
stationary responses (stochastic responses where the statistical mean and variance do not 
change with time) of catchments is an important step to better understanding and 
predicting the often different responses of catchments and their water yields to climate 
(Jones et al., 2012). Individual statistical analyses have emerged that enable us to resolve 
this complexity, including wavelet analysis where decomposition of a time series into a 
time-frequency space enables not only identification of the dominant wavelets (i.e., 
periodicities) in the time series, but also estimation of how the dominant wavelets change 
over time (Torrence and Compo, 1998; Santos and De Morais, 2008). Wavelet cross 
coherence analysis between the spectra of large-scale climatic oscillations and water 
yields can be used to identify specific climatic oscillations responsible for stationary 
signals in the water yield time series. Guiding principles are needed to combine these 
statistical analyses to discriminate non-stationary from stationary signals in catchment 
responses to climate. 

Small, headwater catchments have great potential to serve as sentinels for climate 
change because of their ability to respond rapidly to changing environmental conditions 
including freeze-thaw and flood-drought transitions and trends (Strand et al., 2008). 
Headwater catchments typically have high topographic positions, collect precipitation 
from a relatively small area, and have thin soil layers for moisture storage, making their 
hydrologic regime more sensitive to climate-related changes in environmental conditions 
compared to the higher order catchments into which they drain (Feminella, 1996; Winter, 
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2000; Eimers et al., 2004). Furthermore, headwater catchments contribute a large amount 
of water (Kirby, 1978; Haycock et al., 1993) and cover a large spatial extent (Meyer and 
Wallace, 2001; Sidle et al., 2000; Bishop et al., 2008) of larger catchments, and so have 
considerable influence on what may happen to higher order catchment water yields. In 
addition, these headwater catchments serve as sources of sediments, nutrients, and biota 
to downstream reaches of larger catchment systems (Rabeni and Wallace, 1998; Gomi et 
al., 2002; Wipfli and Gregovich, 2002; Clarke et al., 2008). Therefore, scientific 
investigations of the effects of climate-related environmental changes on water, sediment 
and nutrient export are best made in headwater catchments. 

The purpose of this paper is to present an analytical framework for discriminating 
non-stationary and stationary signals in water yield responses of headwater catchments in 
a forested landscape where climatic variability is evident, but little is known about its 
effect on sustainability of water yields. The (null) hypothesis that climate has no 
significant effect on water yields was tested using a 28-year time series of water yields 
(the longest available within the region) from headwater catchments. If climate does have 
an effect, then a corollary to the hypothesis is that stationary signals are more important 
than non-stationary signals in water yields. Also, if climate does have an effect, then a 
corollary to the hypothesis is that sensitivity to these signals is greater in catchments that 
have lower rates of change in water loading and lower water storage capacity. By 
examining the ability of headwater catchments to respond to both directional climate 
change and non-directional naturally occurring climatic oscillations, it may be possible to 
infer what types of catchments may be the most sensitive sentinels of climate change and 
may see the most extreme impacts of continuing climate warming. 

3.2 Study Area 

Located on the northern edge of the Great Lakes-St. Lawrence Forest Region in the 
Algoma Highlands of Central Ontario (Figure 3.1), the Turkey Lakes Watershed (TLW) 
is about 60 km north of Sault Ste. Marie (47°80’30”N, 84°82’50”W; Figure 3.2). In 
1980 the TLW was established as a long-term experimental research station by Canadian 
federal government agencies to investigate the potential effects of acid rain and climate 
change on terrestrial and aquatic ecosystems (Jeffries et al., 1988). Within the TLW a 
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Figure 3.1: The forest eco-regions of Ontario. 
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Figure 3.2: The Turkey Lakes Watershed. 
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series of headwater catchments representing the range of natural geologic and 
geomorphic variation have been continuously monitored for hydrology and 
biogeochemistry by Natural Resources Canada-Canadian Forest Service and 
Environment Canada scientists. 

The TLW rests on Precambrian silicate greenstone that is metamorphosed basalt, 
with small intrusions of felsic igneous rock (Giblin and Leahy, 1977). The overall relief 
is 400 m, from 644 m above sea level at the summit of Batchawana Mountain to 244 m 
above sea level at the outlet to the Batchawana River. Overlying the bedrock is a thin and 
discontinuous till, ranging in depth from <1 m at higher elevations to 1-2 m at lower 
elevations (Jeffries and Semkin, 1982), although till deposits up to 65 m occasionally 
occur in bedrock depressions (Elliot, 1985). The podzolic soils that have developed in the 
tills follow a generalized sequence of thin and undifferentiated near the ridge, gradually 
thickening, differentiating, and increasing in organic content on topographic benches and 
toward streams (Nicolson, 1988). Highly humified organic deposits occur in wetlands 
(Canada Soil Survey Committee, 1978; Cowell and Wickware, 1983). 

The watershed is covered by a northern tolerant hardwood forest dominated by 
sugar maple ( Acer saccharum Marsh.) (Wickware and Cowell, 1985). There have been 
no disturbances since the 1950s, except an experimental harvest in 1998 on a portion of 
the watershed not included in the present study. Average stand density (904 stems ha" 1 ), 
dominant height (20.5 m), diameter at breast height (15.3 cm) and basal area (25.1 m" ha" 
) are relatively uniform across the uplands, with stand density increasing and dominant 
height decreasing in the wetlands (Jeffries et ah, 1988). White pine ( Pinus strobes L.), 
white spruce ( Picea glauca Moench Voss.), ironwood ( Ostrya virginiana (Mill.) K. 
Koch), and yellow birch ( Betula alleghaniensis Britton) comprise less than 10% of basal 
area in the uplands. The sparse understory of upland stands is dominated (95%) by 
saplings and seedlings of sugar maple as well as a variety of herbs and ferns. Wetland 
stands are mixtures of black ash ( Fraxinus nigra Marsh.), eastern white cedar ( Thuja 
occidentalis L.), red maple (Acer rubrum L.), balsam fir ( Abies balsamea (L.) Mill.), 
yellow birch, and tamarack (Larix laricina (DuRoi) K. Koch.). The understory in the 
wetlands is composed of the seedlings and saplings of overstory trees and various herbs 
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and ferns. There has been no to minimal natural or anthropogenic disturbance to the 
forest cover within these four catchments. 

The watershed is influenced by a continental climate with average 28-year annual 
precipitation of 1200 mm and average annual temperature of 5 °C (Figure 3.3). These 
headwater catchments are influenced by a snowpack that persists from late November, 
early December through to late March, early April and peak flows that occur during 
snowmelt and again during autumn storms. Future climate scenarios for the area predict 
temperatures increasing by 1.0 to 3.0 °C by 2040 and 2.5 to 4.5 °C by 2070, and 
precipitation increasing by 0 to 10% by 2040 and 0 to 20% by 2070 (Price et al., 201 1). 

The four selected catchments, c35, c38, c47 and c50, exhibit a gradient in water 
loading in the form of rain and snow and in water storage in wetlands (percent wetlands 
by area) (Figure 3.4, Table 3.1). To establish the magnitude of the gradient in water 
loading, precipitation records from 1981 to 2010 of six TLW precipitation network 
gauges (namely, the main meteorological recording station outside of the TLW and five 
precipitation gauges distributed within the TLW) were used to calculate the long term 
monthly average precipitation for each of station. The long term monthly average 
precipitations of all the stations were then used to interpolate grid maps of long term 
monthly average precipitation for the whole TLW regions using Anusplin Version 4.37. 
The monthly average precipitation layers were added to derive a grid map of long term 
annual average precipitation for TLW region. For each studied catchments, the long term 
annual average precipitations were calculated by taking the average of all annual average 
precipitation grids that fall within the catchment boundaries. Based on an analysis of the 
precipitation and discharge, the higher elevation catchments had on average 
approximately 5% more precipitation than the lower elevantion catchments (4.16% in c47 
and 5.13 in c50), and on average approximately 20% more discharge (12.34% in c49 
(which was used instead of c47 because of a dysfunctional weir in c47 that resulted in 
significant bypass of water from the weir and therefore an underestimate in discharge, F. 
Beall, Canadian Forest Service, Pers. Comm.) and 26.80% in c50). The additional 15% 
in discharge relative to precipitation may be attributed to underestimates in the 
interpolated precipitation at these specific catchments due to inadequate spatial 
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Figure 3.3: Annual hydroclimatic variation in the Turkey Lakes Watershed region. 
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Figure 3.4: Topographic maps and descriptions of the four catchments selected from the 
Turkey Lakes Watershed for analysis. In c47, there is a buried side-channel, causing loss 
of groundwater (F. Beall, pers. comm). Catchment c49 is a high-elevation catchment that 
has similar average water loading (1327mm/yr) and a relatively small wetland (4%) with 
average discharge of 710 mm/yr, giving a closer approximation to the total discharge 
expected from c47. 
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Table 3.1: Characteristics of catchments in the Turkey Lakes Watershed. 


Characteristics 

c35 

c38 

c47 

c50 

Size (ha) 

4.02 

6.46 

3.43 

9.47 

% Wetland 

1.06 

20.54 

0.36 

10.03 

Elevation of weir (masl) 

386 

415 

503 

507 

Slope (degrees) 

19.36 

13.51 

20.84 

13.45 
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representation of precipitation gauges, or to undetected water inputs from groundwater 
systems passing through fractured bedrock from outside of the catchment boundary. 

To establish the magnitude in the gradient in water storage, we assumed the water 
storage capacity per unit area of the upland soils to be similar. The thin till that overlies 
the bedrock consists of two layer: an upper ablation till that is relatively permeable, and a 
lower basal till which is relatively impermeable. The ablation till that overlies the basal 
till has an average depth of 0.5 m (ranging from <0.2 to 1.0+ m) (Nicolson, 1988). 
However, some of the catchments contain wetlands that may modify water storage 
capacity. For example, catchments c35 and c47 have low water storage potential (1.1% 
and 0.4% wetland, respectively), whereas c38 and c50 contain high water storage 
potential (20.5% and 10.0%, respectively) (Creed et al., 2008). These wetlands likely 
vary in depth, but field measurements indicate that they can have depths up to about 10 m 
(Creed, Unpubl. Data). The modifying affect of these wetlands can either be to increase 
or decrease water export, depending on the antecedent water conditions within the 
wetlands (i.e., when wetlands are full of water, upland drainage waters bypass the storage 
by flowing over the surface of the wetland, and when wetlands are empty of water, 
upland drainage waters are trapped in the wetland as the water table within the wetland is 
below the weir. 

3.3 Methods 

3.3.1 Analytical Framework 

The analytical framework for analysing non- stationary (linear trends) and stationary 
(non-linear oscillations) signals from the computed yearly water yield time series is 
presented in Figure 3.5. 

3.3.2 N on-stationary linear trends 

Annual data were used because daily, monthly and seasonal water yield data were too 
variable and linear trends were not statistically significant (data not shown). Daily water 
yield data from four catchments (c35, c38, c47 and c50) in the TLW were summed to 
derive a yearly time series of water yield running from 1981 to 2008 water years (June- 
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Figure 3.5: Flow chart summarizing the analytical steps. 
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May, with the water year indicated as the calendar year in which the water year began). 
The computed yearly time series of water yields were analysed for non-stationarity, 
which was then removed. The detrended portions were then further studied for the 
presence of multiple stationary signals using wavelet analysis. 


3.3.3 Stationary non-linear oscillations 


Following the method outlined in Torrence and Compo (1998), wavelet power spectra of 
the detrended portion of the yearly water yield data were derived for the catchments 
under investigation. The spectrums were computed by convoluting each time series with 
a scaled and translated version of a transforming wavelet function, usually known as the 
“mother wavelet.” Assuming an equal time interval of At in a time series Xn (n = 0, 1 . . . 

N ~ 1), the corresponding wavelet function, y/ o (rj) (tj is a non-dimensional time 

parameter on which the function depends), must have a zero mean and be localized in 
time and frequency space (Farge, 1992). A Morlet wavelet was chosen due to its 
extensive applications in similar studies involving hydroclimatic time series analysis. In 
addition to its ability to represent the shape of hydrological signals (Kang and Lin, 2007), 
Morlet wavelets provide adequate time and superior frequency resolution compared to 
other wavelet types (Labat, 2005; Soniat et al., 2006). 

The Morlet wavelet is characterized by a Gaussian modulated plane wave: 

y/ 0 (>l) = 71 A e' m ° n e 2 (3.1) 

Where i = an imaginary unit for complex valued Morlet function to overcome 
disadvantages of shift sensitivity, poor directionality and lack of phase information in 
using conventional real valued wavelets; co o is the non-dimensional angular frequency, 
which by default is taken to be 6 to satisfy the admissibility condition (Farge, 1992). The 
continuous wavelet transform of time series X n , for each scale s at all n. with respect to 
the wavelet function y/ 0 (r/) is mathematically represented as (3.2). 
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Where W n ( s ) stands for wavelet transform coefficients, >// for the normalized wavelet, 

(*) for the complex conjugate, s for wavelet scale, n for localized time index, and n! for 
translated time index. 

The transformed signal, W„ (5), is a function of the wavelet scale and translation 
parameters found by conducting the convolution N (the number of data in the time series) 
times for each scale. Other than (3.2), a faster way of performing such a transform is to 
simultaneously calculate all the N convolutions in a Fourier space using discrete Fourier 
transform (DFT) of Xn (3.3). 
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Where k = 0. . . N is the frequency index. In the Fourier space the wavelet transform 
equation (3.2) can be expressed as (3.4). 
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The wavelet power spectrum was calculated as the square of the absolute value of the 
wavelet transform, I W„ (s) I 2 . 

The global wavelet power spectrum (GWPS) based analysis is considered a 
simple and robust technique for characterizing variability in time series (Santos and De 
Morais, 2008). Once the wavelet power spectra were calculated, the GWPS values were 
computed by time- averaging of wavelet spectrum values over all the local spectra. 
Generally, in wavelet based time series analysis, it is believed that the GWPS values 
provide useful information in terms of assessing the scale(s) (or period(s)) contributing 
most to the spectral energy of the time series under investigation. Scales with large 
GWPS values were considered to contribute more spectral energy, while the contribution 
of small GWPS scales were assumed to be little or insignificant. The GWPS values were, 
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therefore, used to identify scale(s) with high potential of periodic signals of stationary 


nature. 


Inverse wavelet transfer of the forward wavelet transform (3.2) above was 
employed to reconstruct signals based on the coefficients of scales identified to have 
dominant GWPS. The inverse wavelet function is mathematically represented as (3.6). 
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Maintaining the coefficients of the chosen scale(s) while replacing all scale coefficients 
by zero and performing backwards (inverse) wavelet transform on the result generated 
the signals characterizing the identified scale. When choosing two or more consecutive 
scales, the inverse transform function builds the reconstructed signal as the sum of 
wavelets of different scales, 5 , at localized time, n. 

Autocorrelation tests were used to examine if there was temporal autocorrelation 
between yearly time series of water yield. The autocorrelation was tested using Matlab’s 
autocorr function and the confidence bounds were found using (3.7). 


nSTDs 
± sqrt (TV) 


(3.7) 


Where nSTDs is the number of standard deviations used in determining significance, and 
N is the number of independent sample data points. Two standard deviations were used to 
approximate 95% significance. 


3.3.4 Searching for multiple stationary signals 

Stationary signals were extracted from the water yield data time series using a step-wise 
methodology involving multiple cycles. The first step of the first cycle used the wavelet 
power spectrum of the detrended water yield time series to identify the dominant 
periodicity in the detrended time series, which was followed by inverse wavelet 
transform to extract the signal of dominant periodicity corresponding to scale(s) of peak 
GWPS values. The extracted signal was then fitted with a periodic function to model and 
statistically analyse the stationarity of the extracted signal. The same steps were followed 
in the second cycle using the residual of the detrended data (detrended water yield minus 
modeled stationary signal from cycle 1) as the input to the wavelet analysis. In the 
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following cycles, subsequent residuals (input data of the previous cycle minus modeled 
stationary signal in the same cycle) were similarily analysed. 

For each cycle, the following “rules” were used to identify existing stationary 
signals: (1) A baseline of GWPS values was determined to select scales with the largest 
GWPS peaks above the set baseline; (2) Peak GWPS values were only considered for 
scales that were able to complete at least two full cycles in the available data (i.e., cycles 
with periods less than 14 years); (3) Subsequent peak GWPS values were analysed until 
their inclusion in the model was not significant; and (4) The same scale could not appear 
twice among signals. After selecting the stationary signals for each catchment, each 
signal was added to the model using forward step-wise regression against the raw water 
yield data, first with the non-stationary signal. 

After independent stationary signals were identified for each nutrient from each 
catchment, an overall stationary model was also developed for each nutrient from each 
catchment by modeling the sum of the independent stationary signals identified. The 
amplitude (difference between maximum and minimum) of the overall stationary signals 
was used to determine the range of change exhibited by stationarity over the study period. 


3.3.5 Wavelet cross coherence analysis 


Wavelet cross coherence analysis was used to understand the dominant oscillations that 
have detectable periodic patterns in the yearly time series of water yield within the 
chosen catchments. Wavelet coherence (Torrence and Webster, 1999; Grinsted et al., 
2004) of two time series X and Y, with wavelet transforms W n x (s) and W„ ¥ (s), is defined 
as 


<G) = 5 - 


(3.8) 


Where S is a smoothing operator both in time and space that can be written as 
S(W) = iS' scalc ( 3'ii me ( W,,(s))) (3.9) 

Where, 5 sca ie denotes smoothing along the wavelet scale axis and 3 lmie for the time axis. 

Equation 3.8 closely resembles that of a traditional correlation coefficient. As a 
result, it is logical to consider wavelet coherence as a localized correlation coefficient in 
time frequency space (Grinsted et al., 2004). The purpose of a smoothing operator is to 
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help locate local maxima or high regions in a plot, such as the cross coherence plot. The 
operation involves a smoothing or bin parameter, which determines the number of nearby 
frequencies, relative to one of interest, averaged (in a manner depending on the 
smoothing filter) to obtain the estimate of the cross coherence at the given frequency, 
producing an estimator that achieves small mean square error, which simultaneously 
keeps the variance and bias small. 

We used wavelet cross coherence to correlate annual average indices for water 
years that were computed from monthly average indices of large scale climatic 
oscillations with the detrended water yield of the chosen catchments at each period using 
a Matlab package from Grinsted et al. (2002). The Multivariate El Nino Southern 
Oscillation (ENSO) Index (MEI), Atlantic Multidecadal Oscillation (AMO), Northern 
Atlantic Oscillation (NAO), Pacific Decadal Oscillation (PDO) indices were selected due 
to their global influence (Figure 3.6, Table 3.2). Wavelet transform plots map the 
correlations of wavelets of varying periods with the water yield time series, and wavelet 
cross coherence plots map the correlations of wavelets of varying periods with two 
different time series, the large scale climatic oscillations and the local water yield time 
series. In both cases, the cone of influence (COI) delimits regions of the plot that can be 
considered as part of the analysis. By analysing the plots, significant associations are 
indicated using warmer colours, and relationships with 95% confidence are circled. For 
wavelet cross coherence plots, arrows indicate the nature of the correlation: arrows to the 
left indicate negative correlations and arrows to the right indicate positive correlations 
(Grinsted et al., 2004). 

3.3.6 Climate drivers of yearly water yield time series 

Correlation tests were performed among the global climate oscillations and between the 
global climate oscillations and local meteorological conditions. To examine correlations 
among the global climate oscillations, Pearson Product Moment correlation tests were 
performed between each of the four climate indices selected. Significant correlations 
were not expected if the global climate oscillations were significantly correlated for a 
short period (i.e., if they were only strongly correlated for 5 years in a 28-year data 
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Figure 3.6: Time series of the Multivariate El Nino Southern Oscillation (ENSO) Index 
(MEI), Atlantic Multidecadal Oscillation (AMO), Northern Atlantic Oscillation (NAO), 
and Pacific Decadal Oscillation (PDO) indices. 
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Table 3.2: Large-scale climatic oscillations used for cross-coherence analysis. 


Indices (Data Source) 

Periodicity 

Climatic influence 

Multivariate El Nino 
Southern Oscillation 
(ENSO) Index (MEI) 
(NO A A - ESRL 
2012a) 

The periodicity of strong ENSO 
events is approximately 2 to 7 
years (Huggett, 1997). 

A positive MEI index (an ENSO 
episode), is associated with warmer 
than normal autumn to spring 
temperatures at the study site 
(Shabbar and Khandekar, 1996, 
Shabbar and Bonsai, 2004). 

Northern Atlantic 
Oscillation (NAO) 
(NCAR 2012) 

There is no statistically 
significant periodicity, although 
there is some evidence of a 7 to 
9 year and approximately 20 
years periodicities (Burroughs, 
2005). 

A positive NAO is associated with 
colder than normal winters at the 
study site (Bonsai et al., 2001). 

Pacific Decadal 
Oscillation (PDO) 
(JISAO 2012) 

PDO tends to occur with a 
periodicity of around 20 to 30 
years (Bum, 2008; Mantua and 
Hare, 2002; Minobe, 1997) 

A positive PDO is associated with 
warmer and drier winters, while 
negative phases are associated with 
cooler and wetter winters at the 
study site (Burn, 2008). 

Atlantic Multidecadal 
Oscillation (AMO) 
(NO A A - ESRL 
2012b) 

AMO is characterised by a 60- to 
90-year oscillation (Knudsen et 
al., 2011). 

A positive AMO anomaly is 
associated with warmer winters and 
lower rainfall at the study site 
(Knight et al, 2006). 
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series). Therefore significant correlations among the global climate oscillations occurred 
only if they existed for nearly the entire time series. 

To examine correlations between the global climate oscillations and local climatic 
conditions, we sought the longest possible meteorlogical records located near the four 
catchments in the TLW: the National Climatic Data Center’s Global Historical 
Climatology Network v2 (GHCN-Monthly) station at Sault Ste. Marie, Michigan (Station 
identifier 42572734000, with a record from 1889 to 2011) and Environment Canada’s 
station at Sault Ste. Marie, Ontario (Station identifier 71260, with a record from 1962 to 
2011). Furthermore, data from the Michigan station were split into two periods (1901- 
1963 and 1964-2010) to compare the relationship between AMO and 
temperature/precipitation during distinct cycles of the estimated 60- to 90- year AMO 
global climate oscillation. 

3.4 Results 

There were significant non- stationary signals in water yields. Catchments showed 
significant negative linear trends in the yearly water yield time series over the past 28 
years (Figure 3.7). The rate of decline in water yield was highest (-14.6 mm/yr) in the 
catchment with lower water loading and lower water storage capacity (c35) (r 2 = 0.52, p 
< 0.001) and nearly as high (-13.0 mm/yr) in the catchment with similar water loading 
but higher water storage capacity (c38) (r 2 = 0.44, p < 0.001). In contrast, the rate of 
decline in water yields was lowest (no significant signal) in the catchment with higher 
water loading and lower water storage capacity (c47) (r 2 = 0.04, p = 0.306), and higher (- 
8.9 mm/yr) in the catchment with higher water loading and higher water storage capacity 
(c50,) (r 2 = 0.18, p < 0.05), but not as high as the lower water loading catchments. Non- 
stationary signals were influenced more by water loading than by water storage capacity. 

There were also significant stationary signals in water yields. The sequence of 
steps for detecting significant stationary signals is illustrated for c35 in Figures 3.8, 3.9 
and 3.10. After removing the linear trend in the water yield data of each catchment 
(Figures 3.8A-3.10A), wavelet transforms (Figures 3.8B-3.10B) and peak analysis 
(Figures 3.8C-3.10C) revealed periods that explained significant amounts of variance in 
the detrended water yield. Sine curves were modeled to the detrended water yield data 
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Figure 3.7: Non- stationary trends in water yield for (A) c35, (B) c47, (C) c38 and (D) c50. 
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Figure 3.8: Stationary trend analysis for the first signal in c35: (A) c35 yearly detrended 
water yield data, (B) wavelet transform (from Equation 3.6) of c35 water yield data 
[strong and weak contributions to the power spectrum are indicated by warm (red) 
colours and cold (blue) colours, respectively, on a logarithmic scale; the thin solid line 
represents the cone of influence, the thick solid lines show the 95% significance level], 
(C) scale selection based on Global Wavelet Power Spectrum [the dotted line at 14 years 
indicates the cutoff for selecting periods of signals (half the 28-year time series)], (D) 
modeling the sine function to the detrended water yield data, and (E) Auto-correlation 
function analysis for the residuals between the observed and modeled water yield. 
Subsequent steps repeated the same process until fitting scales were no longer statistically 
significant. 
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Figure 3.9: Stationary trend analysis for the second signal in c35: (A) c35 yearly 
detrended water yield data, (B) wavelet transform of c35 water yield data, (C) scale 
selection based on Global Wavelet Power Spectrum (D) modeling the sine function to the 
detrended water yield data, and (E) Auto-correlation function analysis for the residuals 
between the observed and modeled water yield. 
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Figure 3.10: Stationary trend analysis for the third signal in c35: (A) yearly detrended 
water yield data, (B) wavelet transform of c35 water yield data, (C) scale selection based 
on Global Wavelet Power Spectrum (D) modeling the sine function to the detrended 
water yield data, and (E) Auto-correlation function analysis for the residuals between the 
observed and modeled water yield. 
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(Figures 3.8D-3.10D), and the autocorrelation function values from the residuals were all 
between 0.5 and -0.5 (Figures 3.8E-3.10E), indicating no significant autocorrelation with 
95% confidence. However, the residuals suggested that periodic characteristics remained, 
so the framework was applied until fitted scales were no longer statistically significant. 

Each of the four catchments showed three significant stationary signals with the 
same periods (i.e., 2 to 2.83 years, 3.36 to 4.76, and 6.73 to 9.51 years), except for c50, 
which had a stationary signal of 5.66 to 8 years instead of 6.73 to 9.51 years (Tables 3.3 
and 3.4). The relative dominance of these three periods varied among the catchments. 

The two catchments with lower water storage capacity (c35, c47) showed the same 
sequence of the 6.73 to 9.51 period detected first, followed by the 2 to 2.83 period, and 
then the 3.36 to 4.76 years period. In contrast, the two catchments with higher water 
storage capacity (c38, c50) showed a different sequence from both the lower water 
storage capacity catchments and each other. 

The 2 to 2.83 year signal was detected first in both c38 and c50, but the sequence 
was 6.73 to 9.51 years, and then 3.36 to 4.76 years in c38 and the opposite in c50 (Table 
3.3). Signals with shorter periods dominated the detrended water yield time series in 
catchments with high water storage capacity, whereas the signals with longer periods 
dominated the time series in the catchments with small wetlands. 

The combination of non-stationary and stationary signals explained the majority 
of variation in water yields for each of the four catchments: 81% (c35), 76% (c38), 56% 
(c47), and 63% (c50) (Figure 3.11, Table 2.4). For catchments with lower water loading 
(c35, c38), the non-stationary signal was stronger than the combination of stationary 
signals. For c35 (low water storage capacity), the non-stationary signal explained 52% of 
the temporal variation in water yields, whereas the stationary signals together explained 
an additional 29%. For c38 (high water storage capacity), the relative importance of the 
two signals shifted slightly, with the non-stationary signal explaining 44% of the 
temporal variation in water yields, whereas the stationary signals together explained an 
additional 32%. In contrast, for catchments with higher water loading (c47, c50), the non- 
stationary signals were lower than the combination of stationary signals. For c47 (low 
water storage capacity), the non-stationary signal explained no variation in water yields, 
whereas the stationary signals explained 56% of the variation in water yields. Similarly, 
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Figure 3.11: Time series for the observed and modelled yearly water yields over the study 
time period for the catchments (A) c35, (B) c47, (C) c38, and (D) c50. 
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Table 3.3: Global wavelet power spectrum (GWPS) values of water yield generated from 
Equation 3.6 in catchments c35 (A), c38 (B), c47 (C) and c50 (D). The detrended water 
yield was used to generate the first signal from the wavelet transform. Water yield used 
for subsequent signals is the difference between water yield from the previous step and 
modeled water yield at the step of investigation. The highlighted values indicate the peak 
GWPS values selected as the stationary signal in each step. See Figures 3. 5-3.7 for 
sample peak selection. 

A - c35 


Level 

Period 

(Years) 

Stationary 
Signal 1 

Stationary 
Signal 2 

Stationary 
Signal 3 

Stationary 
Signal 4 

1 

2 

72,742,838.99 

52,559,388.15 

32,799,491.02 

27,446,013.24 

2 

2.38 

117,046,896.11 

86,142,685.27 

33,099,842.66 

27,794,912.71 

3 

2.83 

95,367,771.17 

71,792,782.62 

6,696,663.58 

4,232,105.67 

4 

3.36 

47,741,733.74 

33,848,471.69 

18,165,496.28 

2,583,822.86 

5 

4 

81,074,024.78 

56,401,908.76 

38,739,832.75 

10,806,888.73 

6 

4.76 

58,234,715.87 

39,741,839.00 

26,980,918.02 

17,925,972.20 

7 

5.66 

51,335,670.69 

20,873,913.33 

14,592,670.07 

10,317,740.44 

8 

6.73 

92,205,150.66 

22,340,644.25 

15,694,400.58 

11,632,192.79 

9 

8 

197,102,724.56 

7,197,672.74 

4,958,396.31 

3,739,473.63 

10 

9.51 

264,112,013.58 

2,417,747.20 

1,526,098.69 

1,211,442.02 

11 

11.31 

45,201,766.78 

9,685,764.83 

6,538,213.21 

5,017,783.08 

12 

13.45 

50,842,798.77 

26,525,867.07 

17,936,846.25 

13,754,007.79 

13 

16 

15,053,551.08 

7,880,583.29 

5,330,643.82 

4,090,358.01 

14 

19.03 

1,053,004.09 

969,869.79 

691,324.14 

586,426.30 

15 

22.63 

9,521,102.60 

9,229,953.56 

6,601,178.77 

5,632,752.84 

16 

26.91 

26,474,669.50 

25,665,094.95 

18,355,442.43 

15,662,607.47 

B - c38 





Level 

Period 

Stationary 

Stationary 

Stationary 

Stationary 


(Years) 

Signal 1 

Signal 2 

Signal 3 

Signal 4 

1 

2 

106,001,656.86 

65,500,872.94 

48,315,829.46 

37,976,716.25 

2 

2.38 

176,776,699.36 

73,227,373.56 

56,593,670.61 

41,443,685.28 

3 

2.83 

163,347,155.50 

29,417,499.54 

24,119,780.39 

11,899,508.67 

4 

3.36 

94,230,821.75 

53,245,045.28 

40,833,464.41 

7,754,080.30 

5 

4 

94,604,196.33 

73,783,452.65 

46,187,212.25 

10,742,035.21 

6 

4.76 

72,194,221.50 

55,545,861.64 

29,425,631.37 

15,231,178.55 

7 

5.66 

63,011,768.80 

51,616,100.99 

17,665,216.18 

9,873,949.32 

8 

6.73 

121,492,134.14 

95,498,143.63 

25,383,305.84 

15,629,556.19 

9 

8 

135,031,783.10 

99,286,457.90 

8,232,465.30 

5,436,313.37 

10 

9.51 

153,895,573.39 

110,475,543.03 

1,763,117.56 

1,617,943.44 

11 

11.31 

18,359,957.45 

12,445,607.56 

2,421,615.86 

1,806,002.61 

12 

13.45 

7,256,840.78 

3,152,322.63 

6,311,266.97 

4,617,283.87 

13 

16 

2,184,458.85 

974,902.91 

1,906,954.17 

1,391,336.79 

14 

19.03 

967,872.18 

978,304.53 

871,358.10 

562,080.93 

15 

22.63 

9,648,040.81 

9,846,927.33 

8,690,551.80 

5,593,546.42 

16 

26.91 

26,827,644.36 

27,380,675.06 

24,165,220.47 

15,553,590.24 


Table 3.3. ( Continued ) 
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C - c47 


Level 

Period 

Stationary 

Stationary 

Stationary 

Stationary 


(Years) 

Signal 1 

Signal 2 

Signal 3 

Signal 4 

1 

2 

157,378,239.90 

130,985,771.78 

103,534,002.71 

64,157,739.24 

2 

2.38 

238,019,826.75 

197,689,749.78 

93,014,746.78 

62,664,581.73 

3 

2.83 

236,765,895.02 

197,011,424.17 

51,846,652.20 

18,944,890.45 

4 

3.36 

193,389,465.55 

160,463,022.85 

105,140,813.67 

10,874,058.26 

5 

4 

197,170,649.88 

162,055,920.48 

129,181,221.50 

4,760,670.53 

6 

4.76 

66,672,825.22 

44,838,963.34 

37,913,995.80 

9,819,489.91 

7 

5.66 

161,946,538.57 

25,973,696.09 

23,481,299.18 

18,331,077.33 

8 

6.73 

423,081,833.54 

45,421,536.92 

39,531,938.84 

28,953,133.67 

9 

8 

207,480,165.06 

81,090,984.09 

62,461,506.46 

46,627,686.57 

10 

9.51 

135,564,307.08 

110,929,919.49 

83,073,858.35 

62,119,444.50 

11 

11.31 

92,175,460.96 

79,159,848.22 

56,239,813.54 

39,268,364.16 

12 

13.45 

220,101,911.49 

190,401,174.75 

134,047,200.73 

92,409,319.83 

13 

16 

65,257,544.17 

56,448,360.73 

39,748,733.61 

27,397,450.38 

14 

19.03 

5,489,763.99 

4,678,540.27 

3,451,728.12 

2,293,780.81 

15 

22.63 

50,655,197.95 

43,105,691.36 

31,948,504.04 

21,155,114.85 

16 

26.91 

140,853,402.31 

119,861,011.50 

88,836,994.23 

58,824,563.32 

D - c50 





Level 

Period (Years) 

Stationary 

Stationary 

Stationary 

Stationary 



Signal 1 

Signal 2 

Signal 3 

Signal 4 

1 

2 

347,530,485.98 

237,768,185.59 

179,026,156.86 

153,497,536.92 

2 

2.38 

392,660,317.07 

129,919,942.14 

91,315,689.24 

78,293,766.94 

3 

2.83 

437,477,169.71 

60,073,959.61 

21,922,598.28 

18,783,311.38 

4 

3.36 

300,915,850.89 

178,837,671.90 

16,154,965.00 

13,633,355.93 

5 

4 

323,010,771.88 

248,900,020.12 

11,055,093.84 

7,544,101.66 

6 

4.76 

110,746,232.35 

78,420,959.73 

18,871,232.68 

7,017,206.14 

7 

5.66 

100,538,905.89 

81,222,487.51 

50,338,414.75 

2,258,900.74 

8 

6.73 

176,988,435.47 

165,444,943.04 

113,468,010.56 

623,426.91 

9 

8 

67,184,433.77 

65,711,576.85 

43,635,232.31 

2,342,173.59 

10 

9.51 

53,122,685.79 

46,737,650.63 

30,147,648.03 

14,541,131.78 

11 

11.31 

309,051,767.56 

215,367,355.62 

153,370,958.01 

122,549,864.32 

12 

13.45 

853,003,199.84 

591,660,421.54 

422,246,007.10 

340,091,344.81 

13 

16 

253,770,934.12 

176,267,394.30 

125,815,125.46 

101,392,846.08 

14 

19.03 

38,213,673.55 

31,454,954.80 

22,839,047.47 

19,529,737.85 

15 

22.63 

368,033,540.57 

305,451,014.96 

221,951,169.09 

190,267,619.61 

16 

26.91 

1,023,365,499.83 

849,346,601.61 

617,164,331.78 

529,064,070.81 
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Table 3.4: Amount of variation explained by and magnitude of change of non-stationary 
(linear trends) and stationary (oscillating cycles) signals in yearly water yield export data 
from catchments c35, c38, c47 and c50 in the Turkey Lakes Watershed. For variation, 
data reported are R" values (with p < 0.05; NS indicating not significant relationships). 
For magnitude, data reported are rate of change in export (mm/year) for non-stationary 
signals and amplitude (mm/year) for stationary cycles. 


Water Export 

C35 

C38 

C47 

C50 

Variation 

Non-stationary signal (%) 

52 

44 

NS 

18 

Trend + Cycle 1 (%) 

65 

59 

19 

39 

Trend + Cycle 1 + 2 (%) 

76 

68 

38 

56 

Trend + Cycle 1 + 2 + 3 (%) 

81 

76 

56 

63 

Stationary signal (%)* 

29 

32 

56 

45 

Magnitude 

Trend (slope) (mm/yr) 

-15 

-13 

NS 

-9 

Cycles (amplitude) (mm/yr) 

±72 

±83 

±79 

±100 
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for c50 (high water storage capacity), the non- stationary signal explained 18% of the 
variation in water yields, while the stationary signals explained an additional 45% of the 
variation in water yields (Table 3.4). The modeled water yields represented the observed 
water yields in each of the four catchments; however, the modeled water yield sometimes 
failed to capture peak water yield, particularly in the catchments with high water loading, 
c47 (e.g., 2003 peak) and c50 (e.g., 1997 peak) (Figure 3.10). 

The wavelet cross coherence plots revealed complex interactions between water 
yield and global climate oscillations (Figures 3.12-3.15). For each of the four catchments, 
there were significant correlations between each of global climate oscillations and water 
yields: (a) MEI was highly correlated to water yields at periods of 2 to 3.5 years, 
particularly in the 1990s (coincident with PDO), with greater MEI resulting in lower 
water yields (Figure 3.12A-3.15A); (b) AMO was correlated to water yields at periods of 
3 to 5 years, although the correlations were stronger at the beginning of the time series 
compared to the end of the time series, with higher AMO resulting in lower water yields 
(Figure 3.12B-3.15B); (c) NAO was correlated to water yields at periods of 3 to 4.5 years, 
particularly in the first half of the time series; however, the nature of the correlation was 
unclear (arrows in different directions) (Figures 3.12C-3.15C); and (d) PDO was highly 
correlated to periods of 2 to 3.5 years, particularly during the five year window from 
1998 to 2003, with higher PDO resulting in lower water yields (Figures 3.12D-3.15D). 

There were significant relationships among major climatic oscillation indices (i.e., 
MEI, AMO, NAO, and PDO) (Table 3.5). The Pearson product moment correlation 
between the annual average values of the climate oscillations revealed significant 
correlations between the MEI and PDO (p<0.001) and the NAO and PDO (p<0.05). Of 
the 12 stationary signals observed in water yield data, six were significantly correlated 
with a climatic index (Table 3.6). MEI was the index most often strongly correlated with 
water yield (seven times total with four correlations being significant). 

Comparing meteorological records at nearby stations with the major climatic 
oscillation indices, only the AMO index was significantly linked to annual temperature 
and precipitation at the TLW meteorological station. Temperature and the AMO index 
were positively linked from 1981 to 2008 (r“ = 0.46, p < 0.001). This positive relationship 
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Figure 3.12: Cross coherence plots using Morlet mother wavelets to generate wavelet 
power spectra (from Equation 1 1) show correlations at each period between yearly water 
yield of c35 and (A) MEI, (B) AMO, (C) NAO, and (D) PDO. The thin solid line 
represents the cone of influence, the thick contours indicate the 95% significance level, 
arrows show the relative phase relationship (left arrows for negative correlations and 
right arrows for positive correlations). Warm colours indicate strong correlations and cool 
colours indicate weak relationships along a linear scale. 
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Figure 3.13: Cross coherence plots using Morlet mother wavelets to generate wavelet 
power spectra (from Equation 11) show correlations at each period between yearly water 
yield of c38 and (A) MEI, (B) AMO, (C) NAO, and (D) PDO. The thin solid line 
represents the cone of influence, the thick contours indicate the 95% significance level, 
arrows show the relative phase relationship (left arrows for negative correlations and 
right arrows for positive correlations). Warm colours indicate strong correlations and cool 
colours indicate weak relationships along a linear scale. 
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Figure 3.14: Cross coherence plots using Morlet mother wavelets to generate wavelet 
power spectra (from Equation 1 1) show correlations at each period between yearly water 
yield of c47 and (A) MEI, (B) AMO, (C) NAO, and (D) PDO. The thin solid line 
represents the cone of influence, the thick contours indicate the 95% significance level, 
arrows show the relative phase relationship (left arrows for negative correlations and 
right arrows for positive correlations). Warm colours indicate strong correlations and cool 
colours indicate weak relationships along a linear scale. 
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Figure 3.15: Cross coherence plots using Morlet mother wavelets to generate wavelet 
power spectra (from Equation 1 1) show correlations at each period between yearly water 
yield of c50 and (A) MEI, (B) AMO, (C) NAO, and (D) PDO. The thin solid line 
represents the cone of influence, the thick contours indicate the 95% significance level, 
arrows show the relative phase relationship (left arrows for negative correlations and 
right arrows for positive correlations). Warm colours indicate strong correlations and cool 
colours indicate weak relationships along a linear scale. 
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Table 3.5: Pearson Product Moment Correlations between the annual average values of 
Multivariant El Nino Southern Oscillation (ENSO) Index (MEI), Atlantic Multidecadal 
Oscillation (AMO), Northern Atlantic Oscillation (NAO), and Pacific Decadal 
Oscillation (PDO) (* indicates p<0.05; *** indicates pcO.OOl). 



MEI 

AMO 

NAO 

PDO 

MEI 

- 

-0.152 

0.011 

0.600*** 

AMO 

- 

- 

-0.387 

-0.314 

NAO 

- 

- 

- 

-0.387* 

PDO 

- 

- 

- 

- 
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Table 3.6: Correlations between climatic indices and the identified stationary signals for 
four catchments in the Turkey Lakes Watershed demonstrating the indices that are the 
most strongly correlated. 


Catchment 

Stationary 

Signal 

Periods 

(years) 

Strongest Correlated 
Index 

r (Direction of 
Correlation) 

P 

c35 

1 

6.73 to 9.51 

PDO 

0.220 (+) 

0.013 


2 

2 to 2.83 

MEI 

0.026 (-) 

0.414 


3 

3.36 to 4.76 

MEI 

0.150 (-) 

0.039 

c38 

1 

2 to 2.83 

NAO 

0.010 (-) 

0.619 


2 

6.73 to 9.51 

AMO 

0.250 (-) 

0.007 


3 

3.36 to 4.76 

MEI 

0.240 (-) 

0.009 

c47 

1 

6.73 to 9.51 

NAO 

0.009 (-) 

0.636 


2 

2 to 2.83 

MEI 

0.014 (-) 

0.552 


3 

3.36 to 4.76 

MEI 

0.210 (-) 

0.015 

c50 

1 

2 to 2.83 

MEI 

0.012 (-) 

0.575 


2 

3.36 to 4.76 

MEI 

0.210 (-) 

0.014 


3 

5.66 to 8 

PDO 

0.054 (-) 

0.232 
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was confirmed using the longer records of the Sault Ste. Marie, Michigan station (1889- 
2010, r 2 = 0.19, p < 0.001) and the Sault Ste. Marie, Ontario station (1962-2010; r 2 = 0.35, 
p < 0.001). Precipitation and the AMO index were negatively linked at the TLW station 
(r = 0.34, p <0.001), though no significant relationships were observed with data from 
the Michigan or Ontario stations. 

Splitting the Michigan meteorological data according to distinct AMO cycles 
(1901-1963 and 1964-2010) revealed a change in the relationships. Split temperature 
records confirmed consistent positive relationships between the AMO index and 
temperature (Figure 3. 16 A), but split precipitation records indicated that from 1901 to 
1963, AMO was significantly positively correlated to precipitation (r 2 = 0.16 ,p < 0.05), 
but AMO was negatively correlated with precipitation during 1964 to 2010 (r = 0.08, p < 
0.05) (Figure 3.16B). 

Looking at the complete AMO cycle from 1901 to 1963 versus temperature 
(Figure 3.16A), the minimum AMO value of -0.39 corresponds to an interpolated 
temperature value of 3.81 °C and the maximum AMO value of 0.36 corresponds to a 
temperature value of 4.96°C. The difference of 1.15°C over 31 years (half the AMO 
period) indicates that 0.037°C per year in climate warming can be explained by a shift 
from a low point in the AMO cycle to a high point, about 55% of the observed 0.067 °C 
per year rate observed in the TLW (Figure 3.3). 

3.5 Discussion 

This paper examined the effects of changing climatic conditions on water yields from 
headwater catchments in a natural (generally undisturbed) northern forested landscape. 
Northern latitudes are expected to see substantial environmental changes based on 
projected climate change scenarios (Carey et al, 2010). We were interested to see if and 
how these effects are already being manifested. By focusing on headwater catchments on 
natural landscapes, we examined the hydrological systems that are most sensitive and that 
have the potential to serve as early warning sentinels of climate change effects to more 
regional catchments, including those that provide water supplies to downstream 
communities. 
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Figure 3.16: Regressions between the AMO Index and the mean annual precipitation (A) 
and temperature values (B) from the Sault Ste. Marie, Michigan weather station during 
two distinct AMO cycles (1901-1964 and 1964-2010). 
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3.5.1 Climate change has effects on water yields 

The analytical framework outlined in this paper (Figure 3.5) enabled the detection of 
climate signals in the water yields of four catchments in the Turkey Lakes Watershed 
(TLW). The framework detected non-stationary linear trends, then after removing these 
trends, detected stationary non-linear cycles caused by climate oscillations. Without 
applying this framework, the non-stationary climate signals would not have been 
identified, discriminated, nor quantified from the original time series of water yields. 

Climate had a substantial non-stationary effect on water yields. A shift to warmer 
and drier conditions resulted in a linear decline in water yields that ranged from 0% (c47) 
to 52% (c35). This decline of yearly water yield may be related to a unidirectional change 
towards warmer and drier conditions over the past 28 years. However, climate also had a 
substantial stationary effect on water yields. Multiple scales of nested periodic cycles 
were detected in the water yield data, explaining variation in addition to the non- 
stationary trends that ranged from c35 (29%) to c47 (56%). The three periodic cycles, 
which occurred from 2 to 2.83 years, 3.36 to 4.76 years, and 6.73 to 9.51 years, were 
found in all catchments, except for c50, where the longest periodic trend was 5.66 to 8.00 
years rather than 6.73 to 9.51 years. A comparison of the periods from the wavelet 
analysis (Figures 3.8-3.10, Table 3.3) and the wavelet cross coherence analysis (Figures 
3.12-3.15) suggests that the 2 to 2.83 year signal was influenced by the MEI and PDO; 
the 3.36 to 4.76 year signal was influenced by the AMO; and the 6.73 to 9.51 year signal 
was influenced partially by the NAO, though it was only strongly correlated to water 
yields for the higher water loading catchments, c35 and c38. A 28-year dataset is 
particularly short for wavelet analysis. It is likely that signals with longer periods are also 
present, but they could not be detected with the data available. 

3.5.2 Climate change results in differential responses among catchments 

There were substantial differences in the non-stationary rate of decline and a marginal 
difference in stationary oscillations in water yield signals among the catchments, despite 
their close proximity (Figure 3.7, Table 3.4). Catchments with lower water loading were 
more sensitive to the non-stationary climate signals (Table 3.4), showing declines in 
water yields of -14.6 mm/year (c35) and -13.0 mm/year (c38). Our findings suggest that 
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there is a higher rate of change in net water loading (precipitation minus 
evapotranspiration and sublimation) on sites that experience less overall water loading, 
leading to a higher rate of change in catchment water yields at low water loading sites 
compared to higher water loading sites. 

Catchments with wetlands formed by topographic depressions and/or flats store 
water (Lindsay et al, 2004; Creed et al., 2003; 2008). Catchments with higher water 
storage capacity buffered the effects of climate warming. The declines in water yield in 
catchments with lower water loading (c35 and c38) were of a greater magnitude than in 
those with higher water loading (c47 and c50). For the catchments with lower water 
loading (c35, c38), a lower water storage capacity increased the amount of variation 
explained in water yields attributed to non- stationary signals (52% for c35 versus 44% 
for c38), but increased the amount of model variation attributed to stationary signals (an 
additional 29% for c35 versus 32% for c38). For the catchments with higher water 
loading (c47, c50), a lower water storage decreased the amount of variation explained in 
water yields attributed to non-stationary signals (0% for c47 versus 18% for c50), and 
increased the amount of model variation attributed to stationary signals (an additional 
56% for c47 and 45% for c50). Higher water storage capacity appeared to reduce the 
effects of shorter periods on discharge, perhaps reflecting the buffering capacity of 
storage capacity on water discharge fluctuations associated with shorter-term periods, but 
not longer term periods. 

Our study indicates that drier (lower water loading) catchments with minimal 
water storage capacity may be more responsive to climate change in humid regions. The 
catchment with the lower water loading and lower water storage capacity (c35) was the 
most sensitive to the non-stationary signal (decline in water yield of -14.6 mm/year) and 
was the catchment for which we were able to explain the most variation (81%) in water 
yields through a combination of the non-stationary and stationary signals. In contrast, the 
catchment with higher water loading and lower water storage capacity (c47), was the 
least sensitive to the non-stationary signal (no significant decline in water yields) and 
only 56% of the variation in water yields could be explained through a combination of 
non-stationary and stationary signals. It is difficult to compare our findings with similar 
studies in the region or elsewhere because there are few studies that investigate the effect 
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of catchment water loading on water yield responses to climate change. Nonetheless, our 
study clearly shows that the climatic characteristics of a region within which a catchment 
is located must be considered prior to selecting the best sentinel for climate change. 

Potential water storage impacted the amplitude of stationary signals, as detailed in 
the conceptual framework in Figure 3.17. In the low water storage catchments (c35 and 
c47), there is minimal storage, so the amplitude of the stationary signals should not 
fluctuate greatly. In the high water storage catchments, low water loading (e.g., c38) 
should result in less discharge because water is stored (bucket empty), while high water 
loading (e.g., c50) should result in more discharge because the wetland is filled (bucket 
full). Accordingly, the amplitude of the stationary signals in high water storage 
catchments is larger, showing a greater response to stationary signals, which could not be 
detected when comparing average water yields (e.g., c35 and c38 in Figure 3.4). 

3.5.3 Global drivers of local responses in water yields 

In the last few decades, scientists have recognized the effects of global climate 
oscillations on regional climate, particularly on the interannual variability of temperature 
and precipitation in regions thousands of kilometres away from the sites of the initial 
fluctuations (Ropelewski and Halpert, 1996; Shabbar el al., 1997; Ionita el al., 2012; 

Zhou et al., 2012). Furthermore, the effects of global climate oscillations on hydrology 
and associated hydroecological processes of water bodies are well documented (Foley et 
al., 2002; Kawahata et al., 2003; Kondrashov et al., 2005; Labat, 2008; Keener et al., 
2010). We contributed to this body of knowledge by exploring the effects of global 
climate oscillations on headwater catchment water yields. 

The non-stationary trend was linked to climate warming. One of the original 
assumptions of this research is that at least part of the observed increase in temperature in 
the TLW (0.067 °C per year, Figure 3.3) is being caused by anthropogenic climate 
warming. We were able to show that the AMO index and temperature are significantly 
correlated both using the TLW data and the longer-term data from the Sault Ste. Marie, 
Michigan weather station over two distinct AMO cycles (1901 to 1963 and 1963 to 2010). 
By using the longer-term data, we were able to show that temperature at the minimum 
AMO values increases by approximately 0.037°C per year until the maximum AMO 
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Figure 3.17: Conceptual model of role of potential water storage on amplitude of 
stationary signals water yield time series. 
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values (approximately 30 years). Coincidentally, the TLW data record (1981 to 2008) 
corresponded with an AMO minimum to maximum half cycle (Figure 3.6), meaning that 
the TLW was being affected both by strong non-stationary and stationary signals. We 
estimate that of the 0.067 °C per year increase in temperature observed from 1981 to 
2008, 0.037°C per year can be explained by the AMO index (from interpolated 
temperatures of 3.8 1°C at the minimum AMO value and 4.96°C at the maximum value 
over a 31-year half-cycle, Figure 3.16A), leaving 0.030°C for other causes, including 
anthropogenic climate change. It is possible that some of the increase in temperature 
ascribed to AMO is also, at least in part, caused by anthropogenic climate change. For 
example, Guan and Nigam (2009) state that non-stationary trends, such as anthropogenic 
climate warming, can affect Atlantic Sub-surface Temperature, potentially aliasing 
effects of global warming into natural climate oscillations. 

The stationary trend was linked to global climate oscillations, with the highest 
frequency (2 to 3.6 years) associated with MEI and PDO, a moderate frequency (3.5 to 6 
years) associated with AMO, and the lowest frequency (6.73 to 9.51) associated with 
NAO. Each of these global climate oscillation indices was correlated to local water yields, 
but the timing, persistence, and direction of correlations varied and only MEI was 
consistent across basins. From the complex patterns in the wavelet cross coherence plots 
(Figures 3.9-3.12), the following generalizations were made: (a) The period of high, 
negative correlations between MEI or PDO and water yields was longer in duration in the 
lower water loading catchments (c35, c38) compared to higher water loading catchments 
(c47, c50) - as the MEI or PDO increase, climate is warmer, drier and water yields 
decrease; (b) The period of high, generally negative correlations between the AMO and 
water yields was shorter in duration in the lower water loading catchments (c35, c38) 
compared to higher water loading catchments (c47, c50) - as the AMO increases, the 
climate gets warmer with less rain and water yields decrease; and (c) The period of high, 
generally positive correlations between the NAO and water yields was longer in duration 
in the lower water storage capacity catchments (c35, c47) compared to higher water 
storage capacity catchments (c38, c50) - as the NAO increases, climate is cooler, wetter 
and water yields increase. Global climate oscillations have prolonged effects on water 
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yields from catchments with lower water loading (c35, c38) and lower water storage 
capacity (c35, c47). 

Some of the global climate oscillation indices had interactive effects. The positive 
correlation between the MEI and PDO indicates that as the MEI increases (reflecting 
warmer, drier conditions) or decreases (reflecting cooler, wetter conditions), the effects 
on climatic conditions may be amplified by the PDO. The negative correlation between 
the NAO and PDO indicates that as the NAO index increases (reflecting cooler, wetter 
conditions) or decreases (reflecting warmer, drier conditions), the effects on climatic 
conditions may be amplified or dampened by the PDO (Table 3.5). 

3.5.4 Future consideration of finer time series of water yield data 

This study used yearly time series of water yields rather than finer time scales. In 
our data analysis, the residuals (observed data minus trend) are consistent with white 
noise, but there may be additional features that this type of data cannot detect, for 
example persistent wetness levels over a period of a few months. Future work using finer 
time scales may reveal if there is a stochastic persistence such as autoregressive time 
series dependence in the noise, or if there is coherence with respect to other climate 
oscillation indices such as NOA or PDO on this finer scale. The cross coherence analysis 
shows there is some correlation or dependence between the catchment data and the global 
climate oscillation indices. In future studies that focus on intra-annual data, we will 
explore if there is lag dependence between the catchment data and previous values of 
these indices. We have not checked the cross coherence of the indices with the residuals 
due to the length of the time series and the fact that the residuals are consistent with white 
noise. This additional step would be interesting to see if the regressions mean estimates 
account for the influence of the climate oscillation indices, or if there is further 
dependence. 

A limitation of the presented analyses (and similar analyses in other studies) is the 
short length of the available data time series. The TLW has one of the longest monitoring 
records in Canada (Creed et al., 201 1); even so, it does not allow for full consideration of 
some of the multi-decadal global climate oscillations. For example, the MEI has a 2- to 7- 
year cycle (Huggett, 1997); NAO has no statistically significant periodicity, although 
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there is some evidence of a 7- to 9-year and approximately a 20-year cycle (Burroughs, 
2005); PDO has a 20- to 30-year cycle (Burn, 2008); and AMO has a 60- to 90-year cycle 
(Knudsen et al. 2011). Complete cycles for MEI and NAO only would be captured in the 
28-year record for the TLW. A potential danger of this short time series is that the 
apparent non-stationary signal may exaggerate the stated influence of climate warming. 
Even still, with an observation window that included only half of the AMO cycle, we 
were able to observe significant correlations between climate and the AMO index and to 
ascribe more than half the non-stationary temperature increase in the TLW to this 
relationship. Continued support and investment in long-term catchment studies, and 
analysis of social factors that influence water supply and demand, will help civic 
institutions adapt to future scenarios of climate and land-use change (Creed et al., 2011; 
Jones et al., 2012). 

3.6 Conclusions 

Climatic change has complex influences on headwater catchment water yields. To 
determine effects of climate change, we must be able to discriminate between 
anthropogenic climate warming (non-stationary trends) and natural climate oscillation 
(stationary cycles) signals on catchment water yields from catchments. Using our 
analytical framework, we were able to show that (A) headwater catchments showed a 
general non-stationary decline in water yields, which may have been caused by the 1 
degree Celsius per decade climate warming that has occurred within the region; however, 
individual catchments varied in their responsiveness to climate change, with catchments 
that have lower water loading and lower water storage capacity being the most sensitive; 
(B) headwater catchments showed stationary cycles in water yields caused by specific 
climate oscillations or possibly the interactive impacts of multiple global climate 
oscillations, which are sometimes significantly correlated with each other; and (C) by 
combining the models of stationary and non-stationary signals identified, it was possible 
to explain the majority of the variance in water yields within each catchment. AMO was 
identified as the index most closely associated with climate in the Turkey Lakes 
Watershed, and it explained approximately 0.037°C of the observed 0.067 °C per year 
increase in temperature from 1981 to 2008. The methodology developed in this study can 
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be used not only to understand hydrologic responses to climate change, but also 
associated biogeochemical responses, including carbon, nitrogen, phosphorus export. 
Future work will extend the analytical framework to include seasonal, monthly and daily 
data as well as predicting both water and solute yields from headwater catchments within 
the TLW and will apply this methodology to other catchments from long-term monitoring 
stations with the longest available records in forested landscapes across northern latitudes 
in North America. 
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Chapter 4: Nutrient export from catchments on 
forested landscapes reveal complex non-stationary and 
stationary climate signals 


Preface 

The previous chapter focused on the role of topography in influencing temporal variation 
in annual water yield, whereas this chapter focuses on the role of topography in 
influencing temporal variation in annual nutrient export. This chapter tests the 
hypotheses that that (1) annual time series of nutrient export contain both non-stationary 
and stationary signals; (2) non-stationary signals are the major signed in terms of 
contributing to high rates of nutrient export over time and may obscure the contributions 
of stationary signals on nutrient export; (3) non-stationary signals of metabolically more 
active nutrients (e.g., nitrogen and phosphorus containing nutrients ) show higher rates of 
change than metabolically less active nutrients (e.g., carbon containing nutrients ); and 
(4) catchments with topographic characteristics that lead to lower hydrologic storage 
potential and lower hydrologic loading potential show the greater rates of changes in the 
trends in nutrient export because the hydrological condition that influence 
biogeochemical processes and nutrient transport are at risk of fundamental changes from 
climate change in these systems. The analytical framework presented in the previous 
chapter is applied to detect trends and cycles in the annual time series of nutrient export 
of carbon (DOC), nitrogen (DON, NOf-N), and phosphorus (TDP) in the same 
catchments. The significance of this chapter is that it uses the analytical framework for 
deconstructing export signeds to revecd changes in trends and cycles in nutrient export 
signeds that will affect the absolute and relative magnitude of nutrient export, which in 
turn will influence downstream aquatic productivity and biodiversity. 

4.1 Introduction 

Scientific evidence for climate change is supported by an increasing body of research 
(IPCC, 2007, Hansen el al., 2012). Climate directly affects water export from catchments 
by controlling the water export (Mengistu el al., 2012), and by extension may also affect 
the export of nutrients dissolved in this water. However, much less is known about this 
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relationship, especially in surface waters (Whitehead et al., 2009). Changes in nutrient 
flow may impact water quality, affecting aquatic productivity and biodiversity, and 
potentially leading to algal blooms that can severely hamper human drinking water and 
recreational activities (e.g., Paerl and Paul, 2012). An understanding of the complex links 
between climate and water resources is needed for planning mitigation and adaptation 
procedures to ensure the sustainability of water resources. 

Climate is inherently dynamic and contains an often complex mixture of non- 
stationary and stationary signals. Non- stationary climate signals have statistical mean and 
standard deviation values that change over time, (e.g., directional climatic warming 
driven by anthropogenic climate change). In contrast, stationary climate signals have 
statistical means and standard deviations that do not change over time (e.g., global 
climate oscillations that are naturally occurring at scales ranging from several years to 
several decades). Discriminating non-stationary signals from the stationary signals that 
characterize natural climatic variability will contribute to a predictive understanding of 
changes in both the quantity and quality of surface waters (Adrian et al., 2009). 

Lakes are a good place to detect changes in climatic conditions because they 
integrate changes that occur within their catchments (Williamson et al., 2008); however, 
headwater catchments may be the best places to detect climate signals for the following 
reasons: (1) they constitute a large spatial extent and contribute a large proportion of 
water to downstream systems (Bishop et al., 2008); (2) they are located at relatively high 
topographic positions, they collect precipitation over a relatively small area, and have 
thin soils for water storage, making their hydrologic regime more sensitive to climate 
related changes in environmental conditions (Strand et al., 2008); (3) their geomorphic 
and geologic complexity influences hydrological flow partitioning (surface/subsurface) 
and pathways tapping into biogeochemical source areas (Creed et al., 1996; Creed and 
Band, 1998); (4) they are found in areas that are usually relatively undisturbed by human 
activities (Blenckner et al., 2007); and (5) there are existing long term monitoring 
networks of headwater catchments particularly in the northern hemisphere (e.g.. Creed et 
al., 201 1; Jones et al., 2012). For these reasons, greater focus should be placed on 
headwater systems and their potential to serve as early warning systems for effects of 
climate change on surface water quantity and quality. 
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In a recent study, we examined catchment water yield responses to climate signals 
in a temperate forest region of eastern North America (Mengistu et al., 2012). We 
developed an approach to remove non-stationary trends and then detect the most 
prominent nested oscillations within the dataset. We found that catchments showed a 
general non-stationary decline in water yields, which may have been caused by the 1°C 
per decade climate warming that has occurred within the region over the past three 
decades; however, individual catchments varied in their responsiveness to changing 
climatic conditions. Catchments with lower water loading and minimal water storage 
capacity were more sensitive to non-stationary and stationary signals both in terms of 
magnitude of change in water yield and proportion of the total variation in year-to-year 
water yield (Mengistu et al., 2012). 

In this study, we extend the analysis presented in Mengistu et al. (2012) to 
identify climate signals within nutrient export records from the same catchments. 

Climate trends and oscillations will affect nutrient export indirectly through alteration of 
water flow partitioning and pathways as well as directly through changes in 
biogeochemistry. This study tested the following hypotheses: (1) yearly nutrient export 
time series for dissolved organic carbon (DOC), dissolved organic nitrogen (DON), 
nitrate-nitrogen (NOf-N) and total dissolved phosphorus (TDP) from TLW headwater 
catchments has significant non-stationary and stationary climatic signals; (2) 
unidirectional non-stationary trends are the major signal as a result of their greater rate of 
contribution to nutrient export over time and may obscure the contributions of stationary 
signals on nutrient export; (3) sensitivity to non-stationary signals is greater for nutrients 
that are more metabolically active (e.g., DON, NOs'-N and TDP) than less metabolically 
active (e.g., DOC) nutrients; and (4) sensitivity to these signals is greatest in catchments 
with relatively low water loading and low water storage capacity, where the hydrological 
system and its hydrological connections that influence biogeochemical processes are at 
risk of fundamental changes from climate change. By examining catchment nutrient 
export responses to changing climatic conditions, it may be possible to infer what types 
of catchments are the most sensitive sentinels of climate change and to evaluate the 
potential magnitude of continued climate warming impacts on dissolved nutrient export 
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to surface waters on forested landscapes. This information could be used to support 
policies of reduced nutrient loads into surface waters. 

4.2 Study Area 

The Turkey Lakes Watershed (TLW) is an experimental forest (47°80’30”N, 
84°82’50”W) located on the Algoma Highlands on the northern edge of the Great 
Lakes-St. Lawrence Forest region (Figure 4.1) and established in 1980 by the Canadian 
federal government (Jeffries et al., 1988). The TLW rests on Precambrian silicate 
greenstone formed from metamorphosed basalt, with small outcrops of felsic igneous 
rock near Batchawana Lake and Little Turkey Lake (Giblin and Leahy, 1977). The 
overall relief ranges from 644 m above sea level at the summit of Batchawana Mountain 
to 244 m above sea level at the outlet to the Batchawana River. A thin and discontinuous 
till overlies the bedrock and ranges in depth from <1 m at higher elevations (with 
infrequent surface exposure of bedrock) to 1-2 m at lower elevations, with till deposits 
up to 65 m occasionally in bedrock depressions. The podzolic soils in the tills are thin 
and undifferentiated near the ridge, gradually thickening, differentiating, and increasing 
in organic content on topographic benches and toward the stream, with highly humified 
organic deposits wetlands (Canada Soil Survey Committee, 1978; Creed et al., 2002). 

The watershed is covered by a northern tolerant hardwood forest dominated by 
sugar maple ( Acer saccharum Marsh.) ( Wickware and Cowell , 1985) with no 
disturbances since the 1950s. Average stand density (904 stems ha' 1 ), dominant height 
(20.5 m), diameter at breast height (15.3 cm) and basal area (25.1 m 2 ha" 1 ) are relatively 
uniform across the uplands, with stand density increasing and dominant height decreasing 
in the wetlands. The uplands overstory includes sugar maple associates, white pine ( Pinus 
strobes L.), white spruce ( Picea glauca Moench Voss.), ironwood ( Ostrya virginiana 
(Mill.) K. Koch), and yellow birch ( Betula alleghaniensis Britton), which comprise less 
than 10% of basal area in the uplands, and the sparse understory is dominated by saplings 
and seedlings of sugar maple and a variety of herbs and ferns. In the wetlands, the 
overstory is a mixture of black ash ( Fraxinus nigra Marsh.), eastern white cedar ( Thuja 
occidentalis L.), red maple (Acer rubrum L.), balsam fir ( Abies balsamea (L.) Mill.), 
yellow birch, and tamarack ( Larix laricina (DuRoi) K. Koch.), and the understory is 
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Figure 4.1: The Turkey Lakes Watershed. 
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composed of the seedlings and saplings of overstory trees and various herbs and ferns 
(Webster et al., 2008). 

The watershed is influenced by a continental climate. The average annual 
precipitation is 1200 mm and the average annual temperature is 5.0 °C [based on a 28- 
year meteorological record (1981-2008) from the Canadian Air and Precipitation 
Monitoring Network (CAPMoN) station located just outside the TLW boundary (see 
Figure 4.1)]. Meteorological conditions within the TLW reflect significant influences 
from Lake Superior in the west and orographic effects created by Batchawana Mountain. 
Snowpack persists from late November, early December through to late March, early 
April, with peaks in stream discharge during snowmelt and again in September to 
November during autumn storms. The average annual atmospheric deposition of sulphur 
and nitrogen [based on a 28-year atmospheric acidic deposition record from the same 
CAPMoN station) show that sulphur deposition has been declining at a rate of almost 1 
kg/ha/yr, but nitrogen deposition has been comparatively stable (Figure 4.2) (M. Shaw, 
Unpublished Data). 

For this study, we examine four catchments, c35, c38, c47 and c50, which exhibit 
a water loading gradient in the form of rain and snow and a water storage gradient in 
wetlands (Figure 4.3). See the conceptual model from Chapter 3 for details on how the 
magnitudes of the gradients in water loading and water storage were determined. In 
summary, for water loading, average precipitation for the four catchmets was interpolated 
from local weather stations in the TLW using Anusplin 4.37. For water storage, we 
compared the wetland size by area of the four catchments and assumed that the water 
storage capacity per unit area of upland soils was similar. 

4.3 Methods 

4.3. 1 Analytical framework 

We applied an analytical framework to identify non- stationary and stationary signals in 
yearly time series of nutrient export from catchments c35, c38, c47 and c50 in the TLW 
(Figure 4.4). Daily, monthly and seasonal nutrient export data were too variable and 
linear trends were not statistically significant (data not shown) so annual data was used 
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Figure 4.2: Nitrogen and sulphur deposition data from 1983 to 2006 (Unpublished data, 
M. Shaw, Environment Canada). 
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4.02 ha 
386 masl 
19.4° slope 

1274 mm/yr average precipitation 
632 mm/yr average discharge 
1.1% wetland 

13.2 kg/ha/yr average DOC export 

1 .0 kg/ha/yr average DON export 

5.2 kg/ha/yr average N0 3 -N export 
0.02 kg/ha/yr average TDP export 



6.46 ha 
415 masl 
13.5° slope 

1257 mm/yr average precipitation 
612 mm/yr average discharge 
20.5% wetland 

12.7 kg/ha/yr average DOC export 
0.8 kg/ha/yr average DON export 

2.7 kg/ha/yr average N0 3 -N export 
0.01 kg/ha/yr average TDP export 
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3.43 ha 
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20.8° slope 

1327 mm/yr average precipitation 
555 mm/yr average discharge 
0.4% wetland 

47.2 kg/ha/yr average DOC export 

2.0 kg/ha/yr average DON export 

1.0 kg/ha/yr average N0 3 -N export 
0.06 kg/ha/yr average TDP export 
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13.5° slope 

1334 mm/yr average precipitation 
776 mm/yr average discharge 
10.0% wetland 

32.5 kg/ha/yr average DOC export 

1.7 kg/ha/yr average DON export 

2.0 kg/ha/yr average NO s -N export 
0.03 kg/ha/yr average TDP export 
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Figure 4.3: Topographic maps and descriptions of the four catchments selected from the 
Turkey Lakes Watershed for analysis. 
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Figure 4.4: Diagram detailing the analytical framework of detecting non-stationary and 
stationary trends in nutrient data. 
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instead. We describe the analytical framework briefly here, but for further details refer to 
Mengistu et al. (2012). 

4.3.2 Yearly nutrient export 

Daily discharge and discharge chemistry data were provided by F. Beall at the Canadian 
Forest Service in Sault Ste. Marie, Ontario. Daily discharges were derived from 
continuously measured stream stage at V-notch weirs on each catchment. The 
concentrations of DOC, NFLfCN, NO 3 -N, and TDP in discharge were derived from 
samples collected every two weeks during the winter, daily during spring snowmelt, and 
weekly or every two weeks during the summer and autumn. Stream water samples were 
filtered through Fisher Q 8 (coarse, fast flow) paper filters; extensive testing found no 
significant difference for samples filtered with fine vs. coarse filters. 

NH 4 + -N, NOa'-N. and total dissolved nitrogen (TDN) were analysed using sodium 
nitroprusside, cadmium reduction and cadmium reduction after autoclave digestion 
methods, respectively, on Technicon Autoanalysers. DON was calculated from TDN 
minus dissolved inorganic nitrogen (DIN = NH 4 + -N + NO 3 -N). DOC was determined by 
removing dissolved inorganic carbon (DIC) by purging with No after acidification, 
converting DOC to DIC by persulfate oxidation catalyzed by UV, then converting the 
resulting DIC to CO 2 by acidification, which was measured by colorimetry. TDP was 
analysed on a Technicon Autoanalyser IIC after autoclave digestion using the 
molybdophosphoric acid colour reaction. 

Daily nutrient concentrations were either measured or estimated by interpolation 
of the two measured values located directly before and directly after the unknown value. 
Daily nutrient exports were then calculated as the product of the total daily discharge and 
the interpolated daily concentrations in the discharge. Yearly nutrient exports were 
calculated as the yearly sum of the daily nutrient exports based on the water year (June 1 
to May 31), before being normalized by catchment area. 

4.3.3 Non- stationary linear trends 

The computed yearly time series of nutrient exports were analysed for non-stationarity by 
linear regression. Significant stationary signals were subtracted from the observed data, 
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and the residuals (i.e., detrended data) were then examined for the presence of stationary 
signals using wavelet analysis. If there was no significance, then no stationary signals 
were removed. 


4.3.4 Stationary non-linear oscillations 


Wavelet transforms display the most prominent frequencies within a yearly nutrient 
export data in a wavelet power spectrum (WPS). WPSs of the detrended portion of the 
yearly nutrient export data were derived using the method outlined in Torrence and 
Compo (1998). A Morlet wavelet was chosen due to its extensive application to similar 
studies involving hydroclimatic time series analysis (Kang and Lin, 2007). 

The WPSs were computed by convoluting each time series with a scaled and 
translated version of a transforming wavelet function. Assuming an equal time interval of 
At in a time series X„(n = 0, 1, . . ., N - 1), the corresponding wavelet function, y/ 0 (>]) 
(where ;/ is a non-dimensional time parameter on which the function depends), must have 
a zero mean and be localized in time and frequency space. The wavelet was characterized 
by a Gaussian modulated plane wave: 

_i _n4_ 

V 0 {n ) = n 2 (4.1) 

Where i = an imaginary unit for complex valued Morlet function to overcome 
disadvantages of shift sensitivity, poor directionality and lack of phase information in 
using conventional real valued wavelets; co 0 is the non-dimensional angular frequency, 
which by default is taken to be 6 to satisfy the admissibility condition (Farge, 1992). 

The continuous wavelet transform of time series X„, for each scale s at all n, with 
respect to the wavelet function y/ 0 (rj) is mathematically represented as 4.2: 


W G) 


— V x ,y/ 


r i 

s 


(4.2) 


Where W„ ( s ) stands for wavelet transform coefficients, t//for the normali z ed wavelet, (*) 
for the complex conjugate, s for wavelet scale, n for localized time index, and ;/' for 
translated time index. 

The transformed signal, W„ (5), is a function of the wavelet scale, and translation 
parameters are found by conducting the convolution N (the number of data in the time 
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series) times for each scale. A faster way of performing such a transform is to 
simultaneously calculate all the N convolutions in a Fourier space using discrete Fourier 
transform of X„ (4.3): 

x k = —Y^x n e-^ nlN (4.3) 

N n = 0 

Where k = 0. . . A is the frequency index. 

In the Fourier space the wavelet transform equation (4.2) can be expressed as 

(4.4): 

W n (s) = Yj x kV * G®* (4.4) 

k = 0 


Where y/(sco) is Fourier transform of a function i//(l / s) and a>k is an angular frequency 
defined as (3.5): 
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The WPSs were calculated as the square of the absolute value of the wavelet transform, 
\W n (s)l 2 . 

Once the WPSs were calculated, the global wavelet power spectrum (GWPS) 
values were computed by time-averaging wavelet spectrum values over all the local 
spectra. GWPS values can be used to assess the scales (or periods) contributing most to 
the spectral energy of the time series under investigation. Scales with large GWPS values 
were considered to contribute more spectral energy, while the contribution of small 
GWPS scales were assumed to be little or insignificant. GWPS values were, therefore, 
used to identify the periods of strongest stationary periodic signals. 

Inverse wavelet transformation of the forward wavelet transform (4.2) was 
employed to reconstruct signals based on the coefficients of scales identified to have 
dominant GWPS. The inverse wavelet function is mathematically represented as (4.6). 
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Maintaining the coefficients of the chosen scales while replacing all scale coefficients 
with zero and performing backwards (inverse) wavelet transform on the result generated 
the signals characterizing the identified scale. For choices of two or more consecutive 
scales, the inverse transform function builds up the reconstructed signal as the sum of 
wavelets of different scales, s, at localized time, n. 

Autocorrelation was used to detect temporal autocorrelation within yearly nutrient 
export data using the autocorr function in Matlab. The confidence bounds were found 
using 4.7: 


nSTDs 
± sqrt ( N ) 


(4.7) 


Where nSTDs is the number of standard deviations used in determining significance, and 
N is the number of independent sample data points. Two standard deviations were used to 
approximate 95% significance. 


4.3.5 Detecting stationary signals 

Stationary signals were detected in the yearly nutrient export data using a step-wise 
methodology. A wavelet power spectrum of the detrended nutrient export data was used 
to identify the dominant periodicity in the detrended time series, which was followed by 
an inverse wavelet transform to extract the signal of dominant periodicity corresponding 
to scales of peak GWPS values. The dominant peak was selected for further analysis only 
if it had not been selected previously, as spectral leaking from one cycle to the next could 
occur. The extracted signal was then fitted with a periodic function to model and 
statistically analyse the stationarity of the extracted signal. The same steps were followed 
in the second cycle using the residuals of the detrended nutrient export data minus 
modeled stationary signal from the first cycle as the input to the wavelet analysis. In each 
of the subsequent cycles, residuals (i.e., input data of the previous cycle minus modeled 
stationary signal in the same cycle) were similarily analysed. 

Subsequent peak GWPS values were considered as long as they made a 
significant contribution to the regression models. After selecting the stationary signals for 
each catchment, each signal was added to the respective nutrient export model using 
forward step-wise regression, starting with the non- stationary signal. After modeling the 


108 



stationary structures at each step for each nutrient, the model at each step was added to 
the overall nutrient model using forward step-wise regression against the raw nutrient 
export data, beginning with the non-stationary component of the nutrient export time 
series. 

An overall stationary model was also developed for each nutrient from each 
catchment by modeling the sum of the stationary structures identified at each step for 
each nutrient. The amplitude (difference between maximum and minimum divided by 
two) of the overall stationary model was used to determine the range of change exhibited 
by stationarity over the study period. 

4.3.6 Wavelet cross coherence analysis 

Wavelet cross coherence analysis was used to determine if there was an association 
between yearly indices of global climate oscillations and the detrended yearly nutrient 
export data from the catchments using a Matlab package based on one from Grinsted et al. 

v 

(2002). Wavelet coherence of two time series X and Y, with wavelet transforms W„ ( s ) 
and W„ Y (s), was defined as 4.8: 

|w„ x G)f)x5G"‘|w;( s 

Where S is a smoothing operator to help locate local maxima both in time. In the plot 
this can be written as 4.9: 

S(W) = S scale (S time (W n (s))) (4.9) 

Where, 5 sca ie denotes smoothing along the wavelet scale axis and 3 lmic smoothing along 
the time axis. 

The Multivariate El Nino Southern Oscillation (ENSO) Index (MEI), Atlantic 
Multidecadal Oscillation (AMO), Northern Atlantic Oscillation (NAO), Pacific Decadal 
Oscillation (PDO) indices were selected due to their global influence. MEI (annualized 
from monthly data based on water year) was used instead of alternative indices [e.g., the 
Southern Oscillation index (Sea Level Pressure differences) and ENSO 3.4 (SST)] 
because it includes six variables: sea level pressure, surface zonal wind, meridional wind, 
cloudiness, sea surface temperature, and surface air temperature (Wolter and Timlin, 


if) 
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1998), which makes MEI effective at discriminating ENSO impacts on hydrology, 
including rainfall and runoff variability (Kiem and Franks, 2001). 

The annual average oscillation indices were calculated from monthly index values. 
Wavelet transform plots mapped the correlations between wavelets of varying periods 
with the yearly nutrient flux time series to readily visualize the dominant frequencies, 
while wavelet cross coherence plots mapped the correlations of two wavelet transformed 
time series (in this case time series of yearly water yield and average annual oscillation 
indices) at local times and scales. The cone of influence delimits regions not influenced 
by edge effect as a result of using time series of finite length, warmer colours represent 
good matching between correlating time series and the thick black solid lines encircles 
regions of the 95% confidence level in relatioship. For wavelet cross coherence plots, 
arrows indicate the nature of the correlation: arrows to the left indicated the two signals 
are in anti-phase, while arrows to the right indicate the two signals are in-phase (Grinsted 
et al., 2004). Angled arrows give an indication of potential lags within the data. Arrows 
pointed downward indicate that X is leading Y by 90°, while upward pointing arrows 
indicate Y leading X by 90° (where 90° represents a time equivalent to one fourth of the 
periodicity represented by the scale of interest). Finally, Pearson Product Moment 
correlation tests were performed between each of the four climate indices studied and 
each of the stationary signal models isolated from within the time series of each nutrient. 
All statistical analyses were completed using Sigmaplot 12.0 (Systat Software, 2012), and 
MatFab version 2010b (Math works, 2010). 

4.4 Results 

4.4.1 Non- stationary and stationary signals 

Significant negative linear trends were detected in DON, NO 3 -N, and TDP export flux 
data, but not DOC export data (Figures 4. 5-4. 8 ). These trends were often strongest in c35, 
the catchment with low water loading and low water storage, with rates of -48 g 
DON/ha/yr, -81 g NOa'-N/ha/yr, and -0.8 g TDP/ha/yr. In contrast, these trends were 
lowest in c47, the catchment with high water loading and low water storage, with rates of 
-22 g DON/ha/yr and no significant trend in NO 3 -N or TDP. However, the weir in c47 
was not ideally placed, and there was a bypass channel that circumvented the weir. The 
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Figure 4.5: Spider plots (A) and line plots (B) showing time series of dissolved organic 
carbon (DOC) export in four catchments in the Turkey Lakes Watershed. 
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Figure 4.6: Spider plots (A) and line plots (B) showing time series of dissolved organic 
nitrogen (DON) export in four catchments in the Turkey Lakes Watershed. 
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Figure 4.7: Spider plots (A) and line plots (B) showing time series of nitrate (NOa’-N) 
export in four catchments in the Turkey Lakes Watershed. 
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Figure 4.8: Spider plots (A) and line plots (B) showing time series of total dissolved 
phosphorus (TDP) export in four catchments in the Turkey Lakes Watershed. 
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closest comparable catchment is c49, which had significant declines in DON and NO 3 -N 
(R 2 = 0.165 and 0.199, respectively) of -21 g DON/ha/yr and -52 g N 03 ~-N/ha/yr, 
respectively. Where nitrogen declines were detected, rates were higher in NO 3 -N than in 
DON. Where DON and TDP declines were detected, rates were higher in DON than TDP 
by almost two orders of magnitude (Tables 4. 1-4.4). 

Significant nested stationary cycles were detected in all nutrient export data and in 
all catchments (Figures 4.9-4.12, Tables 4. 1-4.4). Stationary signals were identified at 
very short (2 to 3 year), short (3 to 5 year), medium (5 to 9 year) and long (9 to 13 year) 
scales (Table 4.5). The very short (very high frequency) 2.0 to 2.8 year signal was 
observed for all nutrients in all catchments. The short (high frequency) 3.4 to 4.8 year 
signal was observed within the DON, NO3 -N and TDP export time series in all 
catchments (except c35 for TDP). The short DOC signal was observed only in c38, the 
low water loading and high water storage catchment. A medium (moderate frequency) 
signal was identified in many catchments, at either 5.7 to 8.0 years or the overlapping 
period of 6.7 to 9.5 years. There was a long signal in all catchments for DOC, and a long 
signal for DON was observed in all catchments except c35. A medium signal was 
observed in only the low water storage catchments (c35, c47) for NO3 -N. There was a 
medium stationary signal for TDP in the low water loading catchments (c35, c38). 
Generally, the earlier cycle (5.7 to 8.0 years) was identified in the nitrogen bearing 
nutrients (i.e., DON, N 03 _ -N)), while the later cycle (6.7 to 9.5 years) was identified in 
the carbon and phosphorus bearing nutrients. The longest cycle identified across the sites 
was 9.5 to 13.5 years and was detected in DON (c50) and in NO3 -N (c35 and c38). 

For each catchment, the overll nutrient export stationary signals revealed 
temporally varying phase relationships (Figures 4.13 and 4.14), at times coupled and at 
time uncoupled, leading to temporal changes in the stoichiometry of nutrient export [i.e., 
temporally varying DOC: nitrogen export ratio and TDP and nitrogen export ratio]. 

4.4.2 Non- stationary vs. stationary signals 

The relative strengths of the non- stationary and stationary signals were evaluated 
first by assessing the magnitude of export attributed to the signal. To assess whether non- 
stationary or stationary patterns contributed more to variability in nutrient export, the 
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Figure 4.9: Wavelet cross coherence plots using Morlet wavelets showing correlations at 
each period between yearly DOC yield and MEI, NAO, PDO, and AMO at four 
catchments in the Turkey Lakes Watershed. The thin solid line represents the cone of 
influence, the thick contours indicate the 95% significance level, arrows show the relative 
phase relationship (left arrows for negative correlations and right arrows for positive 
correlations). Warm colours indicate strong correlations and cool colours indicate weak 
relationships along a linear scale. 
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Figure 4.10: Wavelet cross coherence plots using Morlet wavelets showing correlations at 
each period between yearly DON yield and MEI, NAO, PDO, and AMO at four 
catchments in the Turkey Lakes Watershed. 
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Figure 4.1 1: Wavelet cross coherence plots using Morlet wavelets showing correlations at 
each period between yearly NOa'-N yield and MEI, NAO, PDO, and AMO at four 
catchments in the Turkey Lakes Watershed. 
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Figure 4.12: Wavelet cross coherence plots using Morlet wavelets showing correlations at 
each period between yearly TDP yield and MEI, NAO, PDO, and AMO at four 
catchments in the Turkey Lakes Watershed. 
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Figure 4.13: Final stationary signal models for DOC, and DON export. 
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Figure 4.14: Final stationary signal models for TDP, NOs'-N and DON export. 
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N0 3 -N, DON (kg/ha/yr) N0 3 -N, DON (kg/ha/yr) 


Table 4.1: Amount of variation and magnitude of change explained by non-stationary 
(linear) and stationary (periodic) trends in yearly DOC export data from catchments c35, 
c38, c47 and c50 in the Turkey Lakes Watershed. Data reported are R 2 values (p < 0.05; 
NS indicating not significant relationships), rate of change in export (mm/year) for non- 
stationary trends and amplitude for stationary cycles. 


DOC Export 

C35 

C38 

C47 

C50 

Variation 

Non-stationary signal (%) 

NS 

NS 

NS 

NS 

Trend + Cycle 1 (%) 

29 

34 

20 

18 

Trend + Cycle 1 + 2 (%) 

53 

67 

50 

53 

Trend + Cycle 1 + 2 + 3 (%) 

b 

88 

69 

b 

Stationary Signal (%) a 

53 

88 

69 

53 

Magnitude 

Trend (slope) (g/ha/yr) 

NS 

NS 

NS 

NS 

Trend related overall change (g/ha) 0 

e 

e 

e 

e 

Cycles (amplitude) (g/ha/yr) 

±2,587 

±9,121 

±3,104 

+4211 

Cycle related overall change (g/ha) d 

5,174 

18,242 

6,208 

8,434 


a Amount of variation explained by the final model minus amount of variation explained 
by non-stationary trend. 

b The last stationary signals identified in c35 and c50 DOC export data were not 
significant. 

c Calculated as a product of magnitude of rate of decline (slope) and the number of study 
years. 

d Calculated by multiplying sample amplitude by 2. 

e Trend related overall decline was not calculated as the non-stationary signal identified 
in c38 and c47 NOa’-N export data were not significant. 
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Table 4.2: Amount of variation and magnitude of change explained by non-stationary 
(linear) and stationary (periodic) trends in yearly DON export data from catchments c35, 
c38, c47 and c50 in the Turkey Lakes Watershed. Data reported are R 2 values (p < 0.05; 
NS indicating not significant relationships), rate of change in export (mm/year) for non- 
stationary trends and amplitude for stationary cycles. 


DON Export 

C35 

C38 

C47 

C50 

Variation 

Non-stationary signal (%) 

80 

26 

29 

39 

Trend + Cycle 1 (%) 

86 

43 

50 

51 

Trend + Cycle 1 + 2 (%) 

89 

55 

67 

59 

Trend + Cycle 1 + 2 + 3 (%) 

b 

69 

76 

68 

Trend + Cycle 1 + 2 + 3 + 4 (%) 

b 

b 

b 

73 

Stationary Signal (%) a 

9 

43 

47 

34 

Magnitude 

Trend (slope) (g/ha/yr) 

-48 

-38 

-22 

-43 

Trend related overall change (g/ha) c 

1,248 

988 

572 

1,118 

Cycles (amplitude) (g/ha/yr) 

±137 

±336 

±163 

±221 

Cycle related overall change (g/ha) d 

274 

672 

326 

442 


a Amount of variation explained by the final model minus amount of variation explained 
by non-stationary trend. 

b The last stationary signals identified in c35, c38 and c47 DON export data were not 
significant. 

c Calculated as a product of magnitude of rate of decline (slope) and the number of study 
years. 

d Calculated by multiplying sample amplitude by 2. 

e Trend related overall decline was not calculated as the non-stationary signal identified 
in c38 and c47 NOa’-N export data were not significant. 
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Table 4.3: Amount of variation and magnitude of change explained by non-stationary 
(linear) and stationary (periodic) trends in yearly NO 3 -N export data from catchments 
c35, c38, c47 and c50 in the Turkey Lakes Watershed. Data reported are R” values (p < 
0.05; NS indicating not significant relationships), rate of change in export (mm/year) for 
non-stationary trends and amplitude for stationary cycles. 


NO3 -N Export 

C35 

C38 

C47 

C50 

Variation 

Non-stationary signal (%) 

28 

NS 

NS 

39 

Trend + Cycle 1 (%) 

61 

31 

24 

55 

Trend + Cycle 1 + 2 (%) 

73 

42 

39 

65 

Trend + Cycle 1 + 2 + 3 (%) 

82 

52 

60 

b 

Trend + Cycle 1 + 2 + 3 + 4 (%) 

85 

b 

b 

b 

Stationary Signal (%) a 

57 

52 

60 

26 

Magnitude 

Trend (slope) (g/ha/yr) 

-81 

NS 

NS 

-52 

Trend related overall change (g/ha) c 

2,106 

e 

e 

1,352 

Cycles (amplitude) (g/ha/yr) 

±740 

±213 

±521 

±413 

Cycle related overall change (g/ha) d 

1,480 

426 

1,042 

826 


a Amount of variation explained by the final model minus amount of variation explained 
by non-stationary trend. 

b The last stationary signals identified in c38, c47 and c50 NCR'-N export data were not 
significant. 

L Calculated as a product of magnitude of rate of decline (slope) and the number of study 
years. 

d Calculated by multiplying sample amplitude by 2. 

e Trend related overall decline was not calculated as the non-stationary signal identified 
in c38 and c47 NO 3 -N export data were not significant. 
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Table 4.4: Amount of variation and magnitude of change explained by non-stationary 
(linear) and stationary (periodic) trends in yearly TDP export data from catchments c35, 
c38, c47 and c50 in the Turkey Lakes Watershed. Data reported are R 2 values (p < 0.05; 
NS indicating not significant relationships), rate of change in export (mm/year) for non- 
stationary trends and amplitude for stationary cycles. 


TDP Export 

C35 

C38 

C47 

C50 

Variation 

Non-stationary signal (%) 

42 

16 

NS 

19 

Trend + Cycle 1 (%) 

51 

49 

14 

38 

Trend + Cycle 1 + 2 (%) 

62 

71 

32 

49 

Trend + Cycle 1 + 2 + 3 (%) 

b 

77 

b 

b 

Stationary Signal (%) a 

20 

61 

32 

30 

Magnitude 

Trend (slope) (g/ha/yr) 

-0.8 

-0.9 

NS 

-0.6 

Trend related overall change (g/ha) c 

20.8 

23.4 

e 

15.6 

Cycles (amplitude) (g/ha/yr) 

±4 

±12 

±3 

±5 

Cycle related overall change (g/ha) d 

8 

24 

6 

10 


a Amount of variation explained by the final model minus amount of variation explained 
by non-stationary trend. 

b The last stationary signals identified in c35, c47 and c50 TDP export data were not 
significant. 

c Calculated as a product of the magnitude of rate of decline (slope) and the number of 
study years. 

d Calculated by multiplying sample amplitude by 2. 

e Trend related overall decline was not calculated as the non-stationary signal identified 
in c38 and c47 NCL’-N export data were not significant. 
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Table 4.5: Observed stationary signals for DOC, DON, NO3 -N, and TDP in four 
catchments in the Turkey Lakes Watershed. 


Nutrient 

Signal 

Periods for observed stationary signals in each catchment (years) 

c35 

c38 

c47 

c50 

DOC 

1 

6.7 -9.5 

6.7 -9.5 

6.7 -9.5 

6.7 -9.5 


2 

2 . 0 - 2.8 

2 . 0 - 2.8 

2 . 0 - 2.8 

2 . 0 - 2.8 


3 


3.4 -4.8 

4.0 -5.7 


DON 

1 

2 . 0 - 2.8 

5.7 -8.0 

6.7 -9.5 

9.5-13.5 


2 

3.4 -4.8 

2 . 0 - 2.8 

2 . 0 - 2.8 

5.7 -8.0 


3 


3.4 -4.8 

3.4 -4.8 

2 . 0 - 2.8 


4 




3.4 -4.8 

NO 3 --N 

1 

5.7 -8.0 

2 . 0 - 2.8 

3.4 -4.8 

3.4 -4.8 


2 

9.5-13.5 

3.4 -4.8 

5.7 -8.0 

2 . 0 - 2.8 


3 

3.4 -4.8 

9.5-13.5 

2 . 0 - 2.8 



4 

2 . 0 - 2.8 




TDP 

1 

6.7 -9.5 

6.7 -9.5 

2 . 0 - 2.8 

3.4 -4.8 


2 

2 . 0 - 2.8 

2 . 0 - 2.8 

3.4 -4.8 

2 . 0 - 2.8 


3 


3.4 -4.8 
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magnitude of declines (maximum minus minimum) in nutrient export associated with the 
linear trend was compared with the range between amplitudes of the nutrient export 
associated with the overll stationary signal models (Tables 4. 1-4.4). For DOC, no 
significant trends were identified at any catchments and therefore the non-stationary 
signal was less strong than the stationary signals and dominated variability within the 
DOC export time series. For DON, NO.f-N, and TDP, the magnitude of the overall 
change in the non-stationary signal was larger than the range of variation between the 
amplitudes of the overall stationary signal (Tables 4. 2-4.4). For DON, all catchments had 
non-stationary signals that accounted for more of the range in the observed export fluxes 
than any of the identified stationary frequencies. Stationary cycles resulted in greater 
export of DON from the catchments with high water storage capacity (c38, c50) than 
those with low water storage (c35, 47). For NOa"-N, where trends existed, the magnitude 
of change in nutrient export explained by non-stationary signals was higher than that 
explained by overall stationary signals. For TDP, signals were more complex. The 
magnitude of change explained by c38 non-stationary and overall stationary models was 
almost the same, and there was no significant non-stationary signal for c47, while the 
magnitude of the non-stationary signal change in c35 and c50 was larger than the overll 
stationary signal models. 

The strength of the non-stationary and stationary signals was also assessed by 
examining the magnitude of variation in export attributed to each signal (Tables 4. 1-4.4). 
In most catchments, stationary signals explained more variation in nutrient export than 
non-stationary signals. The exceptions were c35 (low water loading, low water storage), 
where 0.80 of the total 0.89 variation in DON (i.e., 90%) and 0.42 of the total 0.62 
variance in TDP (68%) and c50 (high water loading, high water storage), where 0.39 of 
the total 0.73 variation in DON (53%) and 0.39 of the total 0.65 variation (60%) in NO3- 
N, were due to the non-stationary signal. For all other catchments and nutrients, the non- 
stationary signals explained less of the total explained variation in nutrient export (Table 
4. 1-4.4). Looking at the nutrients, 53% to 0.88% of the variation in DOC concentrations 
was explained by stationary signals and no significant non-stationary signal was 
identified in any of the catchments. The other nutrients behaved more similarly to each 
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other. For DON, 0.09 to 0.47 of the variation was explained by stationary signals; for 
NO3 -N, 0.20 to 0.61 of the variation was explained by stationary signals; and for TDP, 
0.26 to 0.60 of the variation was explained by stationary signals. Therefore, other 
nutrients had weaker, but still important stationary signals. 

4.4.3 Combined signal strength 

Coefficient of determination values were derived by correlating the linear trend model 
added to the sine wave oscillatory models with the raw data. The combined non- 
stationary and stationary models were able to explain 66% of variation within DOC, 77% 
for DON, 66% for NO3 -N, and 55% for TDP (Tables 4. 1-4.4) when averaged across the 
four catchments. In general, the combined stationary and non- stationary signals were able 
to explain the most variation in nutrient export in catchments with low water loading 
(0.72 for both c35 and c38). 

4.4.4 Linking climate and nutrient signals 

Several stationary signals were correlated with global oscillation indices (Table 4.6). 
Several correlations were not significantly correlated at p<0.05, but are presented to show 
the general relationships among the catchments and the global oscillations. For DOC, the 
most prominent stationary signals were most often correlated with the AMO, though only 
one catchment was significantly correlated. For nitrogen species, the predominant 
stationary signals were most often correlated with the NAO. For TDP, the results are 
more heterogeneous and complex than those for the other nutrients, with no clear climate 
oscillation found to be strongest across catchments. 

The very short 2.0 to 2.8 year and the long 9.5 to 13.5 year signals were not 
significantly correlated to any of the tested global climate oscillations (Table 4.6). The 
short 3.4 to 4.8 year signals were correlated with the MEI (DON and TDP) and NAO 
(DON and NO3 -N) indices. The medium 5.7 to 8.0 and 6.7 to 9.5 year signals were 
correlated with AMO (DOC and TDP), NAO (DON and N0 3 '-) and PDO (TDP). 
Generally, the earlier cycle (5.7 to 8.0 years) was identified in the nitrogen bearing 
nutrients (i.e., DON, N03 _ -N)), while the later cycle (6.7 to 9.5 years) was identified in 
the carbon and phosphorus bearing nutrients. 
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Table 4.6: Results of Pearson correlation tests between isolated stationary signals (short = 3.4 to 4.8 year cycle; medium = 5.7 to 8.0 
or 6.7 to 9.5 year cycle) and climate oscillation indices. Several correlations were not significantly correlated at p<0.05, but are 
presented to show the general relationships among the catchments and the global oscillations. 


Signal 

c35 

Index r 2 (Dir) P 

c38 

Index r 2 (Dir) P 

c47 

Index r 2 (Dir) P 

c50 

Index r 2 (Dir) P 

DOC 

Short 

Medium 

AMO 0.20 (-) 0.065 

AMO 0.28 (-) 0.024 


AMO 0.20 (-) 0.065 

DON 

Short 

Medium 

ME1 0.16 (-) 0.042 

MEI 0.20 (-) 0.024 

NAO 0.13 (-) 0.073 

NAO 0.22 (+) 0.016 

NAO 0.15 (-) 0.049 

NAO 0.18 (+) 0.033 

NOf-N 

Short 

Medium 

NAO 0.22 (-) 0.015 

NAO 0.14 (+) 0.056 

NAO 0.15 (+) 0.052 

NAO 0.14 (+) 0.057 

TDP 

Short 

Medium 

PDO 0.14 (+) 0.063 

MEI 0.12 (-) 0.089 

AMO 0.22 (-) 0.015 

MEI 0.19 (-) 0.025 

MEI 0.17 (-) 0.038 


129 




4.4.5 Differential sensitivity of catchments 

Catchments did not respond in a homogeneous manner to environmental changes 
occurring within the region. The two catchments with lower water loading showed 
greater variability in nutrient export signals. For example, c35 (low water loading, low 
water storage) generally had the strongest linear trends [i.e., largest magnitude of linear 
decline in nutrient export flux and largest coefficient of determination values (Table 4.1- 
4.4)] for nutrients that showed trends. In contrast, c38 (low water loading, high water 
storage) generally had strongest stationary signals for dissolved organic nutrients [i.e., 
stationary signals with largest nutrient export amplitudes and largest coefficient of 
determination value (Table 4. 1-4.4)], but the amplitude of the stationary signal was 
among the smallest for NO 3 -N even though c38 NOs'-N stationary signals has modest 
coefficient of determination to explain variability in NO.f-N export. A combination of 
catchments is needed to capture climate effects on nutrient export - c35 was the most 
sensitive to directional climate change, while c38 was most sensitive for detecting 
climatic oscillations (Table 4. 1-4.4). 

4.5 Discussion 

The effects of climate change on water quality and quantity are not clearly understood 
(Whitehead et al., 2009). Contributions of non- stationary and stationary climate signals to 
nutrient export were discriminated in four catchments on a forested landscape. There are 
few studies that focus on the links between nutrient export from headwater catchments 
and that discriminate between climate trends vs. multiple climate oscillations, so in this 
study we highlight the patterns, and suggest processes to explain them (including those 
that influence nutrient export through changes in hydrological flow partitioning and 
pathways and/or biogeochemical rates of reactions), which can be tested in future studies. 

4.5.1 Climate signals in nutrient export 

Both non- stationary and stationary signals were detected in nutrient export in the Turkey 
Lakes Watershed over a 26-year period. Non-stationary signals were declines in nutrient 
export (with the exception of DOC which had no directional change), while (multiple) 
stationary signals had the effect of increasing or decreasing the magnitude of the decline 
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in nutrient export. Specifically, the combination of signals led to larger declines in 
nutrient export during the warmer/drier portions of the oscillations (reflected in positive 
MEI, positive AMO and/or negative NAO), and to smaller declines in export during the 
cooler/wetter portions of the oscillations. There was one exception, c35, which had low 
water loading, low water storage, and where warmer, drier conditions lead to further 
increases in some of the nutrients, NOf-N and TDP. This finding is supported by 
previous studies that show that during cooler, wetter water years, hydrological flow paths 
may be expanded into new biogeochemical source areas (e.g., vertical (up the soil 
column) and horizontal (up the contributing hillslopes into the outer extents of the 
variable source area)) flushing nutrients from near surface and surface soils, while in 
warmer, drier water years, hydrological flow paths may be contracted (e.g., Creed et al., 
1996; Creed and Band, 1998). The consequence is that global climate warming appears 
to be leading to reduced nutrient export that can be further reduced when oscillations 
enter a warm, dry period. 

The multiple stationary signals occurred at very short (2.0 to 2.8 yr), short (3.4 to 
4.8 yr), medium (either 5.7 to 8.0 yr or the overlapping period of 6.7 to 9.5 yr), and/or 
occasionally at long (9.5 to 13.5 yr) frequencies. The very short signal was observed for 
all nutrients in all catchments but was not significantly correlated to any of the tested 
global climate oscillations as the periodicity of these signals was very short compared to 
the periodicity of the tested climatic oscillations. The short signal was correlated with 
both the MEI and NAO. This signal falls within the observed MEI periodicity (Huggett, 
1997), and the NAO does not have a statistically significant periodicity (Burroughs, 
2005). MEI and NAO were also found to be significantly correlated with each other 
(Mengistu et al., 2012), which may suggest a common climate oscillation influencing the 
catchments. The medium signal was identified in many catchments, at either 5.7 to 8.0 
years or the overlapping period of 6.7 to 9.5 years and was most often correlated with the 
AMO. The periodicity of the AMO is approximately 60 to 90 years (Knudsen et al., 
2011), but previous work has shown that the AMO is strongly correlated with 
temperature in the Algoma Highlands and may explain about half of the observed 
temperature increase in the area over the past three decades (Mengistu et al., 2012). The 
long signal was not significantly correlated to any of the four global climate oscillations. 
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These findings suggest that multiple oscillations are controlling nutrient export and that 
some are linked to global climatic oscillations and others are not (and remain 
unexplained) in part because longer time series are needed. 

4.5.2 N on-stationary signals obscure the contributions of stationary signals 

Compared to stationary signals, non- stationary signals contribute more to the overall 
change in nutrient export dynamics due to having higher rates of decline in nutrient 
export resulting greater cumulative changes over time. For catchments with observed 
non-stationary signals, the magnitudes of decline associated with the non- stationary 
signals during the study period were compared with the ranges between the amplitudes of 
the overall stationary signals (Tables 4.1- 4.4). Non-stationary signals obscure the 
contributions of stationary signals as non-stationary signals revealed greater over all 
decline in the magnitude of nutrient export (approximately 572 to 1,248 g/ha for DON; 
1,352 to 2,106 g/ha for N03~-N; 16 to 23 g/ha for TDP) compared to the ranges explained 
by the amplitudes of the overall stationary signals (approximately 274 to 672 g/ha for 
DON: 426 to 1,480 g/ha for N0 3 "-N; 6 to 24 g/ha for TDP). 

4.5.3 Differential sensitivity of metabolically labile vs. recalcitrant nutrients 

Both non-stationary and stationary climate signals were observed to affect nutrient export 
fluxes. First, for non-stationary signals, Mengistu et al. (2012) reported an average 
decline in water yields of 9.25 mm/yr among these same catchments, but no significant 
negative linear trends were detected in DOC export flux data (Table 4.1). The fact that 
water yield declined but DOC in discharge did not suggests that climate has a 
concentrating effect on DOC export: less water is being exported over time but the 
amount of DOC being exported is staying the same or even increasing relative to the 
amount of water being exported. This is supported by reports of trends for increasing 
DOC concentrations in streams draining forested headwater catchments both in North 
America (Eimers et al., 2008) and Europe (Worrall and Burt, 2007; Sarkkola et al., 2009). 
These studies suggest that changes in stream flow have an important role on trends in 
DOC concentration. In contrast to DOC, significant negative trends were detected in 
DON, NO3 -N, and TDP export flux data of most catchments (except c38 and c47 N0 3 "-N, 
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and c47 TDP) (Tables 4.2- 4.4). Dillon and Molot (2005) observed a much greater 
decrease in total phosphorus and DON compared to DOC in Dorset, Ontario, so 
mobilization of these nutrients from wetlands is sensitive to changes in water table 
position and redox levels. It is possible that DON, NOi'-N and TDP are more labile than 
DOC (Stepanauskas et al., 2000; Kaushal and Lewis, 2005; Lpnborg and Alvarez - 
Salgado, 2012; Hendrickson et al., 2002), and that with warmer, drier conditions, these 
more labile nutrients are more preferentially metabolized instead of being exported. If 
DON, NO3 -N and TDP are preferentially metabolized, then this finding suggests that 
there is a shift from more to less labile DOM entering the surface waters as a result of 
climate warming trends (i.e., shift towards higher C/N ratios). 

Second, the overall magnitude of declines in nutrient export associated with the 
linear trend was compared with the magnitude of the range between the amplitudes of the 
nutrient export associated with the overall stationary signal model (Tables 4.1- 4.4). For 
DOC, no significant trends were identified in any catchments and therefore the non- 
stationary signal was less strong than the overall stationary signal that dominated 
variability within the DOC nutrient export, which ranged from approximately 5, 174 g/ha 
to 18,242 g/ha (Table 4.1). For catchments with observed non- stationary signals, the 
magnitudes of decline associated with the non- stationary signals (approximately 572 to 
1,248 g/ha for DON; 1,352 to 2,106 g/ha for NO3VN; 16 to 23 g/ha for TDP) were larger 
than the ranges explained by the amplitudes of the overall stationary signals 
(approximately 274 to 672 g/ha for DON: 426 to 1,480 g/ha for NO3 -N; 6 to 24 g/ha for 
TDP) (Tables 4.2- 4.4). If non- stationary signals are linked to global warming [even if 
only 50% at most, as was the case for water yield reported in Mengistu et al. (2012)], 
then this suggests that anthropogenic global warming trends have had a greater role in 
reducing labile (N, P) nutrient loads than natural global warming oscillations at multiple 
time scales of years to decades. Generally, catchments with high water storage potential 
(c38 and c50) appeared to show larger amplitudes in the overall stationary signal model 
of DOC, DON and TDP exports compared to catchments of lower hydrologic storage 
potential. However, the overll stationary signal from these catchments revealed smaller 
amplitude for NO;/-N export. 
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In general, nitrogen and phosphorus species appear to be more responsive (more 
labile) than carbon to both climate trends and climate oscillations. This may be due to 
several reasons: (1) there are greater concentrations of nitrogen and phosphorus in surface 
and near-surface layers (shoots and roots, litter fall) where climate reactions dominate at 
the land-surface interface (Devito et al., 2000; Liu et al., 2003; Marland el al., 2003; 
Macrae et al., 2005); (2) there could be higher rates of transformation at this land-surface 
interface due to warmer temperatures; (3) there could be more redox switching from 
oxidizing to reducing environments due to rising water tables and/or expanding/ 
contracting surface flows (Creed et al., 201 1); and/or (4) nitrogen and phosphorus 
compounds may have greater mobility along the surface/near-surface hydrologically 
connected flow pathways (Creed et al., 1996; Devito et al., 2000). These changes have 
implications for the stoichiometry of exported nutrients and therefore on the productivity 
and biodiversity of downstream waters (Seitzinger et al., 2005; Justic et al., 2005; 

Donner and Scavia, 2007). Exported nutrient stochoimetry is also influenced by 
oscillation linked catchment responses in nutrient transport that results temporally 
dynamic phase relationships leading to a variable ratio between ecologically important 
nutrients parameters (Figures 4.13, 4.14). Catchments with lower water loading capacity 
appeared have extreme stoichiometry, c35 and c38 catchments respectively revealing 
lowest and highest stoichiometry in both C:N and P:N ratios. Future studies will focus 
not only on DOM but also on DOM constituents to investigate if there are differences 
among labile vs. recalcitrant DOC pools. Measuring DOC properties (e.g., SUVA, 
Fluorescence Index) might reveal a biogeochemical signature that would provide insight 
into fate of the DOC (e.g., Thurman, 1985; McKnight et al., 2003; Cory et al., 201 1). 

4.5.4 Differential sensitivity among the catchments 

Catchments did not respond in a homogeneous manner to environmental changes 
occurring within the region. Despite a small range in annual average water loading (up to 
10%), climate had a greater effect on nutrient export signals in catchments with lower 
water loading. The catchment with low water loading and low water storage capacity, c35, 
generally had the largest declines in nutrient export flux and the most variation explained 
by non-stationary signals for nutrients that showed trends (Tables 4. 2-4.4). The 
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catchment with low water loading and high water storage, c38, generally had stationary 
signals with the largest nutrient export amplitudes and the most variation explained, 
although the magnitude of stationary signal variations were among the smallest for NO 3 - 
N export (Tables 4. 2-4.4). 

Water storage capacity has a potentially strong effect on nutrient export by 
influencing water residence time in the catchment during which biogeochemical 
transformations can occur that may alter the fate (i.e., storage vs. atmospheric or aquatic 
export) of nutrients. For example, organic matter decomposes and is metabolized quickly 
during the snowfree season in catchments with low water storage capacity 
(Ponnamperuma, 1972; Sahrawat, 2004; Moore et al., 2008). In contrast, organic matter 
can be trapped in catchments with high water storage capacity and decomposition slows 
down. When the water table drops, oxidation may lead to preferential removal of more 
labile nutrients, and when the water table rises, the remaining more recalcitrant nutrients 
will become exported as storage capacity fills and water flows through the catchment 
(Webster et al., 2008; Creed and Beall, 2009; Creed et al., 2012). 

4.5.5 Catchments as sentinels of changing nutrient exports 

Are headwater catchments good sentinels of climate change? To answer this 
question, we consider both the predictability and vulnerability of nutrient export in 
headwater catchments under changing environmental conditions. For predictability, in 
most cases, the non- stationary signals explained less variation in nutrient export than 
stationary signals, suggesting that global climate warming had a smaller (lesser) influence 
on nutrient export than global climate oscillations. However, in c35, the catchment with 
low water storage and low water loading, the non-stationary signal explained more 
variation, suggesting that catchments with few to no wetlands that receive relatively less 
water may be most sensitive to climate change. On this humid landscape, catchments 
with less water loading and storage appeared to be more sensitive to changes in climatic 
conditions. By examining the total variance of the nutrient export explained by stationary 
and non-stationary signals, it is possible to understand what catchments were most 
sensitive to climate. Combined signals tended to explain more of the nutrient export 
variance in catchments with lower water loading potential, a similar pattern to that 
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observed in water yield (Mengistu et al., 2012). Combined signals for nitrogenous 
nutrients were strongest in c35, while c38 had the highest amount of variance explained 
for DOC and TDP. For vulnerability, some catchments showed substantial variation in 
nutrient export over the observed period. Catchment c35 showed the greatest range in 
inorganic nutrient export, NO f-N. with stationary signals varying by as much as 1480 g 
NCV-N/ha (average annual export for this catchment is 3500 g NCff-N/ha/yr, Creed and 
Beall , 2009) and trends showing declines of about 81 g NC> 3 ~-N/ha/yr (for a total of more 
than 2000 g NCb'-N/ha over the observed period). This catchment has little accumulation 
of organic matter because of the steep slopes, shallow soils and minimal water storage 
capacity. In contrast, c38 showed the greatest range in the export of organic nutrients, 
DOC, DON and TDP. Minimum and maximum values of the stationary signals ranged 
by as much as 18,242 g DOC/ha, 672 g DON/ha and 24 g TDP/ha, and linear trends 
showed declines of 40 g/ha/yr for DON and 0.9 g/ha/yr for TDP (no significant trends for 
DOC). 

Comparing these results with observed stationary and non- stationary signals in 
water yield (Mengistu et al., 2012), the non-stationary signal explained more variation in 
magnitude than the overlll stationary signals in c35 for nutrients with observed significant 
non-stationary signal, showing that there is a link between water yield and nutrient export, 
but there are also differences. In nearly every case, for catchments where no non- 
stationary trend was detected, the change in nutrient export related to stationarity was 
much greater than the variation explained by the non-stationary trends. This suggests that 
surface water quality is already shifting away from historical ranges controlled by natural 
fluctuations, and this could be very important for ecosystems adjusted for the levels of 
the past. 

To effectively monitor the effects of both non-stationary and stationary signals of 
climate on headwater catchments, a combination of catchments is needed to capture 
climate effects on nutrient export. The catchment with low water loading and low water 
storage, c35, was the most sensitive to directional climate change, while the catchment 
with low water loading and high water storage, c38, was most sensitive for detecting 
climatic oscillations. This suggests that catchments with relatively low water loading and 
low water storage capacity are good sentinels for detecting the effects of anthropogenic 
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climate change on nutrient export. In contrast, catchments with relatively low water 
loading but high water storage are good sentinels for detecting natural climatic 
oscillations and their effect on nutrient export. 

4.6 Conclusions 

Climate warming and changes in atmospheric circulation that leads to climatic 
oscillations have resulted in clear reductions in water and nutrient export from the Turkey 
Lakes Watershed over the past three decades. Complex links were revealed between non- 
stationary and stationary changes in climatic conditions versus nutrient export to surface 
waters. More variation was explained in labile nutrients (DON, NO 3 -N, TDP), and these 
nutrients were also more sensitive to climate signals. Overall, the catchment with low 
water storage potential and low water loading potential was the most sensitive to non- 
stationary and stationary climatic conditions, suggesting that it was the most effective 
sentinel for a dynamic climate. Although we used one of the longest data sets available in 
Canada, even longer sets would allow the detection of signals with longer periodicities. 
Applying the framework to other catchments in the Turkey Lakes Watershed and around 
the world will determine the broad implications of these results. 
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Chapter 5: Conclusion 


5.1 Research Findings 

This dissertation used a combination of spatial analysis, statistical modeling and wavelet 
analysis to investigate the relationship between topography and the spatial and temporal 
controls of water and nutrient export from small headwater catchments in the Turkey 
Lakes Watershed (TLW) that are dominated by otherwise uniform vegetation and 
bedrock geology. The major contributions of this thesis are the (1) understanding of the 
role of various topography related catchment spatial characteristics on the export of 
nutrients; (2) explanation of the mechanism of topography regulated hydrologic control 
on the export behaviour dissolved organic carbon (DOC), dissolved organic nitrogen 
(DON), nitrate -nitrogen (NO3 -N), and total dissolved phosphorus (TDP); (3) 
development of an analytical framework to detect non-stationary (linear trends) and 
nested stationary (non-linear oscillations) signals from catchment yearly water yield and 
nutrient export responses; (4) understanding of the relative sensitivity of TLW headwater 
catchments to climatic forcing by quantifying the variability explained by non-stationary 
and nested stationary components of catchment responses; and (5) identification of the 
dominant climatic oscillations that have detectable periodic patterns in the yearly time 
series of water yield and nutrient export parameters. 

Chapter 2 tested the hypothesis that topography was the most organizing template 
for hydrology influenced nutrient export control from thirteen relatively undisturbed 
headwater catchments in the TLW. Several topography related metrics of hydrological 
flushing potential (e.g. VSA, effV SA) and hydrological storage potential (e.g. flats, 
depressions) were derived following the methods outlined in Creed and Beall (2009), 
which were then statistically assessed for their role in explaining the observed spatial 
variance in nutrient export behaviours among the studied catchments. Topographic 
metrics representing hydrological flushing and hydrological storage explained most of the 
spatial variance in nutrient export. Higher redox potential (oxygenation) in effVSAs when 
water tables are relatively low favours nitrification (oxidation of NH3 + -H to NOf-N), 
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making effVSA a net source of NCV-N. The strong positive association observed 
between effVSA and catchment NO 3 -N export suggests flushing of accumulated NO 3 -N 
from the effVSA as the main mechanism of NCV-N export when the water table rises and 
intersects the soil layer. EffVSA also showed moderately strong but negative association 
with export of organic forms of nutrients (DOC, DON) and TDP, suggesting continued 
mineralization of organic matter and adsorption of organic matter by mineral soils before 
transported to receiving water bodies. 

Metrics reflecting hydrological storage, however, showed significant positive 
relationships with DOC, DON and TDP, suggesting that an increase in the relative 
proportion of depression and flats in a catchment would increase the export of organic 
forms of nutrients due to presence of hydrologic features favouring mechanisms of 
organic matter production in the catchment. There was a significant negative relationship 
between metrics reflecting hydrologic storage and NO 3 -N export, suggesting that 
wetlands positioned along the hydrologic pathway can intercept NCV-N loaded runoff 
before reaching streams, becoming sinks of NCV-N by facilitating its conversion to 
organic matter or its alternative route in the form of gas efflux from denitrification. 

The focus of Chapter 2 was to investigate the relationship between topography 
and temporal responses to climate signals in water yield time series from the TLW. The 
sensitivity of four topographically and hydrologically representative headwater 
catchments were analyzed to test the hypothesis that annual water yield time series are 
composed of a complex mixture of non-stationary linear and cyclic stationary climatic 
trends that are manifested differently in catchments with different topographic 
organization. An analytical framework was developed to reveal if the selected catchments 
had differential sensitivity to non-stationary climate change and naturally occurring 
stationary climate oscillations in water yield time series over 28 years. The variability 
explained by non-stationary signals from annual water yield time series were 
independently analyzed, and wavelet analysis was employed to disassociate the 
remaining stationary component into signals of various periodicities. 

The analytical framework detected non-stationary (linear) declines in all 
catchment water yield data except catchment c47 and three major stationary cycles in the 
water yields in all the studied catchments. The observed non-stationary catchment 
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responses may have been linked to the unidirectional global warming of one-degree 
Celcius per decade rise in temperature in the region. Stationary responses might be linked 
to climate oscillations, including AMO, NAO, PDO and ENSO. Each of these global 
climatic oscillation indices showed correlation with catchment water yield, even though 
the timing, persistence, and direction of correlation varied. Despite the limitation of 
having only 28 years of water yield data, cross coherence analysis revealed that AMO 
showed the strongest influence in the regional hydrology, with consistently high negative 
correlations. AMO also exhibited strong correlations with TLW climatic variables, 
suggesting it might partially explain the unidirectional global warming experienced in 
TLW region. 

Generally, catchments with different topographic organizations exhibited 
differential sensitivity to climate change signals. Catchments with lower water loading 
and lower water storage potential showed a stronger decline in water yield over time 
compared to catchments with higher water loading and higher water storage potentials. 
The catchment with the lowest water loading and the lowest water storage potential (c35) 
was the most sensitive to non- stationary signals, while the catchment with the highest 
water loading and the lowest water storage potential (c47) was the least sensitive to non- 
stationary climatic signal, suggesting that non-stationarity is affected more by water 
loading potential than water storage potential. Catchments with higher water storage 
potential, however, appeared to buffer fluctuations in water yield associated with shorter- 
term periods. 

The main focus of Chapter 4 was investigating topographic influences on climate 
related catchment temporal responses in nutrient export. The analytical framework that 
was developed and tested in Chapter 3 was applied to detect non- stationary and stationary 
signals in yearly time series of DOC, DON, NO3 -N, and TDP nutrient export. There was 
a significant negative non-stationary (linear) trend in all nutrients except DOC. The lack 
of significant declines in DOC coupled with the significant declines in water yield 
declines suggest that there is a concentrating effect of climate on DOC. Significant nested 
stationary signals were also detected in all nutrient export data and in all catchments at 
very short (2 to 3 year), short (3 to 5 year), medium (5 to 9 year) and long (9 to 13 year) 
scales. Some of the observed stationary signals were linked to global climatic 
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oscillations, but others were not, in part because longer time series are needed to detect 
oscillations of longer periodicities. 

In general, the combined stationary and non-stationary signals were able to 
explain most of the variation in nutrient export, with more variation explained in 
catchments with low water loading potential. The non-stationary linear trend explained 
the most variation in the catchment with relatively low water loading and low water 
storage capacity, suggesting these traits are characteristic of good sentinels for detecting 
the effects of anthropogenic climate change on nutrient export. In contrast, stationary 
trends explained the most variation in the catchment with relatively lower water loading 
but higher water storage potential, suggesting these catchment topographic organizations 
may be good sentinels for detecting natural climatic oscillations and their effect on 
nutrient export. More variation was explained in metabolically more active or labile 
nutrients (DON, NO^'-N, TDP) than metabolically less active nutrients (DOC), and these 
nutrients were also more sensitive to climate signals. 

While wavelet and wavelet coherence analyses are becoming more widespread 
when studying hydrological time series and relationships between time series, this 
chapter takes an interesting and slightly different approach than previous wavelet-based 
investigations in hydrology by de-convoluting dominant signals in the time series in an 
attempt to understand the predominant climate forcing. Without applying this framework, 
these climate signals would not have been identified, discriminated, nor quantified from 
the original water yield and nutrient time series. 

5.2 Recommendations for Future Research 

Landscape physical characteristics play important roles in controlling the natural 
hydrological and biogeochemical dynamics in catchments (Neff et al., 2003; Christopher 
et al., 2008; Deng et al., 2011), and therefore not only control spatial differences in 
catchment nutrient export but also influence climate-linked catchment temporal responses 
in water yield and nutrient export. This thesis examined the role of topography in 
controlling spatial and climate linked temporal export responses of water and nutrient 
from relatively pristine headwaters. The following suggestions will continue to improve 
our understanding of these complex systems: 
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1. Conduct similar investigations in geographic areas of different climate and/or 
long term monitoring stations: A 28-year data set is particularly short for time 
series and wavelet analysis. Trends in the data may be missed or ascribed to non- 
stationary signals. Longer time series will allow the detection of signals with 
longer periodicities. 

2. Extend the analytical framework to investigate finer time scales (seasonal, 
monthly and daily data sets): This will allow us to see if signals of stationary 
responses reveal themselves when short frequency time series are used for both 
water yield and nutrient export investigations. Multiple stationary signals can be 
detected in shorter frequency time series, even if non- stationary linear trends 
could be difficult to detect because of the inherent seasonal variability in finer 
scale time series. While annual time scale investigations are important to 
understand yearly nutrient dynamics and long term trends, annual investigation 
may mask finer scale biogeochemical processes strongly influencing the overall 
biogeochemical dynamics in catchments (Oczkowski et al., 2006). 

3. Explore the effectiveness of the tested topographic metrics in explaining nutrient 
export at larger scales and in different forested regions: As a catchment scale 
increases so does the heterogeneity in landscape spatial characteristics, which in 
turn has implications for associated catchment responses, making scale the most 
important issue in process based investigations. Different processes become 
important at different spatial scales and processes that are important at small 
scales may not be important at large scales as new dominant processes may 
emerge as the scale gets larger (Meentemeyer and Box, 1987; Bloschl and 
Sivapalan, 1995). One way of understanding process behavior in large scale 
catchments is to investigate small-scale process behavior in hopes of extending 
that behavior to larger scales (Frisbee et al., 2012). Catchments also vary in their 
biogeochemical responses based on forest type (Amiotte-suchet et al., 2007), 
forest disturbance (Lamontagne et al., 2000; Smith et al., 2012) and deposition 
regime (Sleutel et al., 2009). Extending this investigation to other catchments of 
varying scale, physical, environmental and landuse conditions is necessary in 
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order to have a broader understanding of the controls on catchment nutrient 
export. 

4. Explore the effectiveness of the tested topographic metrics in predicting nutrient 
export at finer (seasonal, monthly and weekly) time scales: Such exploration 
would confirm if the tested metrics are capable of capturing finer time scale 
processes, for example seasonally distinct biogeochemical processes, which are 
important to larger scale export phenomena. This would further our process based 
understanding of nutrient export controls. 

5. Extend the spatial analysis to evaluate the effectiveness of the tested topographic 
metrics in predicting other forms of nutrients at various time scales: The extension 
would enhance our understanding of the underlying linkages between nutrients 
and the mechanisms of their export control. This would, in turn, have implications 
for the management of watershed for water protection. 

5.3 Anticipated Significance 

The research developed a new framework for analysing the role of climate (unidirectional 
climate change and stationary climatic variations) in catchment responses. The 
framework employed a novel approach by de-convoluting dominant signals in catchment 
responses in an attempt to understand the predominant climate forcing. Without applying 
the proposed framework, these climate signals would not have been identified, 
discriminated, nor quantified from the original time series of water yields and nutrient 
export responses. The framework can, therefore, be extended to other environmental 
studies involving time series analysis. Furthermore, the study introduced new spatial 
metrics that are robust in terms of modeling nutirent export fluxes from the studied 
catchments in addition to testing the spatial metrics outlined in Creed and Beall (2009) on 
a range of nutrients having longer time series. 

This research improves the scientific understanding of catchment (hydrological 
and biogeochemical) responses to various temporal and spatial environmental factors. 
There is considerable climatic variability and significant unidirectional climate change in 
the Algoma Highlands region, and this study provided a key insight in revealing the 
complex links between directional climate change, large scale oscillations of atmospheric 


146 



circulation, and water yield and nutrient fluxes from a range of topographically, and thus 
hydrologically, diverse but relatively pristine forested catchments. In addition, the study 
increased our understanding of the most important topography-related catchment 
characteristics on the export of range of different nutrients of organic and inorganic 
nature. Furthermore, the investigation provides comprehensive information about natural 
nutrient dynamics and export behaviors against which similar conditions in disturbed 
situations can be compared. 

The Algoma Highlands region is characterized by granite-dominated bedrock 
geology that makes water bodies and aquatic ecosystems in the region susceptible to the 
effects of acid rain as a result of low buffering capacity (Scheider et al., 1979; 
Shrivastava, 2007). Without well-planned land-use management strategies, climate 
change and increased rates of environmental disturbance in the region from resource 
extraction threatens not only the quality and quantity but also the ecosystem integrity of 
water resources in the region. Practically, the knowledge gained from the research can be 
used to support policies of sustainable development in the region at present and in the 
future. 
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Appendix 1: Evaluatating the role of topography in 
nutrient export using two windows of investigation 


The appendix investigates the role of topography in explaining the observed spatial 
variation in nutrient export among the headwater catchments of Turkey Lake Watershed 
(TLW) using two independent periods of investigation. The first window of investigation 
run from water year 1981 to 1996 and involved all the 13 headwater catchments. The 
second window of investigation extended the study period to 2008 water year, but 
removed all the catchments that experienced disturbance in 1997 forest harvesting 
experiment (catchment c31, c33 and c34). Topographic metrics discussed in Chapter 2 
were regressed against long-term annual average nutrient flux for the two windows of 
investigation. This analysis confirms the results obtained in chapter 2, which used the 
annual nutrient fluxes (kg/ha/a) for the water year from 1981 to 2008 (except for the three 
harvested catchments, c31, c33 and c34, for the period 1981 to 1996 water years). 

Regression analyses examining the relationship between topographic metrics and 
annualized N, P and C export are presented in Table A.l and Table A. 2 for 1981 to 1996 
and 1981 to 2008 windows of investigation respectively. The investigation confirms that 
metrics representing flushing mechanism and flat areas are major players in controlling 
nutrient export and each of these metrics have completely opposite reactions towards 
NO3 -N and DOC, DON or TDP nutrients. NO3 -N showed a strong positive correlation 
with metrics of flushing mechanism while a general pattern of direct association was 
observed between wetlands and DOC, TDP and DOC nutrients. As indicated in Table 
A. 3, most catchment variances in annualized DOC, DON and TDP export were best 
explained by regression models created using W:effVSA as independent variable. 
Annualized NO;f-N exports, however, were well predicted by metrics representing 
flushing mechanism (as indicated by the multiple regression model involving NWeffVSA 
and d 2 effVSA / dTI n 2 ). 
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Table Al.l: Relationships between topography related metrics and annualized nutrient export (1981 to 1996 water years) of the 13 
headwater catchments in the Turkey Lakes Watershed (Direction = Directions of relationship, R 2 = coefficient of determination, and 
NS = not significant, all other relationships p < 0.05). 


Topographic Metrics 

NO 3 -N 

D< 

TC 

DON 

TDP 


Direction 

R 2 

Direction 

R 2 

Direction 

R 2 

Direction 

R 2 

Wetlands (%) 

- 

0.674 

+ 

0.916 

+ 

0.640 

+ 

0.769 

VSA (%) 

NS 

— 

NS 

— 

NS 

— 

NS 

— 

efTVSA (%) 

+ 

0.518 

- 

0.456 

- 

0.526 

- 

0.471 

NWeffVSA (%) 

+ 

0.738 

- 

0.741 

- 

0.679 

- 

0.648 

d-c//VSA/dTI n 2 

+ 


NS 

— 

NS 

— 

NS 

— 

W:VSA 

- 

0.637 

+ 

0.935 

+ 

0.709 

+ 

0.848 

W:effVSA 

- 

0.626 

+ 

0.937 

+ 

0.729 

+ 

0.869 

NWeffVSA and d 2 c//VS A/dTI n 2 

+ 

0.873 

— 

— 

— 

— 

— 

— 
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Table A1.2: Relationships between topography related metrics and annualized nutrient export (1981 to 2008 water years) of the 13 
headwater catchments in the Turkey Lakes Watershed (r = correlation, R 2 = coefficient of determination, and NS = not significant, all 
other relationships p < 0.05). 


Topographic Metrics 

NO 3 -N 

D< 

TC 

DON 

TDP 


Direction 

R 2 

Direction 

R 2 

Direction 

R 2 

Direction 

R 2 

Wetlands (%) 

- 

0.630 

+ 

0.920 

+ 

0.766 

+ 

0.814 

VSA (%) 

NS 

— 

NS 

— 

NS 

— 

NS 

— 

efTVSA (%) 

+ 

0.710 

- 

0.496 

- 

0.552 

- 

0.497 

NWeffVSA (%) 

+ 

0.855 

- 

0.815 

- 

0.779 

- 

0.706 

d 2 e//VSA/dTI n 2 

+ 


NS 

— 

NS 

— 

NS 

— 

W:VSA 

- 

0.632 

+ 

0.930 

+ 

0.819 

+ 

0.888 

W:effVSA 

- 

0.638 

+ 

0.934 

+ 

0.841 

+ 

0.909 

NWeffVSA and d 2 e//VS A/dTI n 2 

+ 

0.887 

— 

— 

— 

— 

— 

— 
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Table A1.3: Best regression models for predicting average annual C, N and P export from TLW headwater catchments. Two models, 
one for each window of investigation, are presented for each parameter. 


Response Variable 

R 2 ,p 

Equation 

N0 3 "-N (kg N ha 1 a' 1 ) 

0.873, <0.001 

N0 3 -N = 0.966 + (0.241 x NWeffW SA) + (49969.886 x d 2 e#VSA/dTI 2 n ) 

N0 3 "-N (kg N ha 1 a 1 ) * 

0.887, <0.001 

N0 3 -N = 0.892 + (0.274 x NWeffVSA) - (25521.513 x d 2 e#VSA/dTI 2 n ) 

DOC (kg C ha 1 a 1 ) 

0.937, <0.001 

DOC = 13.519 + (10.541 x W: effW SA) 

DOC (kg C ha 1 a 1 ) * 

0.934, <0.001 

DOC = 13.052 + (10.527 x W: effW SA) 

DON (kg N ha 1 a' 1 ) 

0.729, <0.001 

DON = 1.113 + (0.328 xW: effW SA) 

DON (kg N ha 1 a' 1 ) * 

0.841, <0.001 

DON = 0.895 + (0.352 x W: effWSA) 

TDP (kg P ha 1 a 1 ) 

0.869, <0.001 

TDP = 0.0158 + ( 0.0123 x W: effVSA) 

TDP (kg P ha 1 a 1 ) * 

0.909, <0.001 

TDP = 0.0143 + ( 0.01 18 x W: effW SA) 


* Equations for 1981 to 2008 window of investigation. 
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be in the form of a check or money order referencing your account number and this 
invoice number RLNK500921072. 

Once you receive your invoice for this order, you may pay your invoice by credit 
card. Please follow instructions provided at that time. 

Make Payment To: 

Copyright Clearance Center 
Dept 001 
P.O. Box 843006 
Boston, MA 02284-3006 

For suggestions or comments regarding this order, contact RightsLink Customer 
Support: customercare @ copyright.com or +1-877-622-5543 (toll free in the US) or 
+1-978-646-2777. 

Gratis licenses (referencing $0 in the Total field) are free. Please retain this 
printable license for your reference. No payment is required. 
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